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COMPUTER ANALYSIS 


OF RING-STIFFENED SHELLS OF REVOLUTION 
By Gerald A. Cohen 

Structures Research Associates, Laguna Beach, California 


SUMMARY 


This report presents the equations and method of solution for a 
series of five compatible computer programs for structural analysis of 
axisymmetric shell structures. User manuals and other program 
documentation for these programs are presented in a separate companion 
report. These programs, designated as the SRA programs, apply to a 
common structural model but analyze different modes of structural response. 
They are: 

(1) Linear asymmetric static response (SRA 100) 

(2) Buckling of linearized asymmetric equilibrium states (SRA 101) 

(3) Nonlinear axisymmetric static response (SRA 200) 

(4) Buckling of nonlinear axisymmetric equilibrium states (SRA 201) 

(5) Vibrations about nonlinear axisymmetric equilibrium states (SRA 300) 

The theory of a sixth related program, for the imperfection sensitivity 
analysis of buckling modes of nonlinear axisymmetric equilibrium states, 
has been presented in a previous NASA report. 

The structural model treated is a branched shell of revolution with 
an arbitrary arrangement of a large number of open branches but with at 
most one closed branch. The shell wall is assumed to be of orthotropic 
material with principal axes of orthotropy in meridional and circumfer- 
ential directions. Geometric properties of the structure may vary only 
in the meridional direction; material properties of the shell wall may 
vary in the thickness direction as well as the meridional direction. Also 
treated are: 

(1) discrete irotropic ring attachments, 

(2) isotropic stringers, whose stiffness is circumferentially 
distributed, and 

an idealized elastic foundation attached to the shell wall. 


( 3 ) 



INTRODUCTION 


During the past decade an almost bewildering variety of computer 
programs has been developed for the analysis of shell structures (ref . 1) . 
When one narrows the field to those designed for elastic shells of 
revolution, he is still confronted with the names of at least forty 
authors in this country alone who have been active in developing pro- 
grams of overlapping capabilities (refs. 1 and 2). At the time reference 
2 was written, however, there were known to be only four major systems 
which cover the most common problems of stress, buckling, and vibration 
of elastic shells of revolution. In addition to the SRA programs of 
this report, these include two finite-difference programs, BOSOR (ref. 3) 
and SALORS (refs. 4 and 5), and Kalnins’ forward Integration programs 
(ref. 6). 

The SRA programs employ the Zarghamee version of the forward 
integration method (ref. 7) for the solution of the basic linear bound- 
ary value problem. This method requires the calculation of only four 
complementary solutions, as opposed to the usual eight, over open branches. 
The main features of the present system of programs which have not been 
generally available in the other systems are: 

(1) buckling analysis under general asymmetric loads, 

(2) imperfection sensitivity analysis, and 

(3) branched shell capability (see fig. 1) 

User documentation for the present system of six programs is 
presented in a companion report (ref. 8). As these programs have been 
developed over a period of time, the theory underlying some of them has 
already been published in the open literature (refs. 9-12). The theory 
of the nonlinear axisymmetric response program and the buckling program 
for general asymmetric equilibrium states, which is a new program, have 
not been previously presented. The purpose of this report is to bring 
together the underlying equations and (improved) method of solution for 
each of these programs except the imperfection sensitivity program, the 
theory of which has been presented in a previous NASA report (ref. 13). 
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SYMBOLS 


r (i) (i) (i) 

W »^2 » c 12 


e 1> e 2 , e i2 


e ,e 
x’ y 


e 

z 


© 1 » e 2» e 12 


F ,F ,F. 

* y ♦ 


e 1»F2» e 3» e 4 


FLi ,FL2 *1X13 , FL4 


GJ 

G«> 

I 


X' 



I 

xy 


K 


^ 1 >^2 »^3 

Li ,L2 


M ,M ,M. 

* y ♦ 


Mi ,M 2 ,Mi2 


ring or stringer cross-sectional area 
ring centroidal radius 

shell wall normal stiffnesses, eqs. (A- 3 ) 

ring or stringer elastic modulus 

orthotropic elastic moduli 

ring centroidal eccentricities relative to 
corresponding boundary point on shell reference 
surface 

normal eccentricity of stringer centroid 
relative to shell reference surface 

linearized shell stretching strains 

effective ring force loads per unit of 
circumferential length 

equivalent shell forces, eqs ( 78 ) 

equivalent ring forces, eqs. ( 79 ) 

ring or stringer torsional stiffness 

shell wall shear stiffnesses, eqs. (A- 3 ) 

stringer section moment of inertia about 
circumferential centroidal axis 

ring section moments of inertia 

structural stiffness 

elastic foundation moduli 

effective shell moment loads per unit ©f area 
ring stress couples 
shell stress couples 
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(i) 


N ,N ,N 
x’ y’ 


n 

P,Q,S 


P 

Rl ,R2 
r 

s,4>,z 


Ti,T 2 ,Ti2 

T <i> 

U 


u,v,w 


w x ,w y ,W <j> 
Xi,X 2 ,X 3 


x,y 

/N 

z 


e l> e 2» e 12 

E <J» 


mas s coefficients for shell inertial loads, eq. (109) 
number of stringers 

effective ring moment loads per unit of 
circumferential length 

circumferential harmonic number 

effective shell forces per unit of circumferential 
length in axial, radial, and circumferential 
directions, respectively 

local pressure for live pressure field at unit X 
meridional and circumferential radii of curvature 
small circle radius 

meridional, circumferential, and normal 
coordinates, respectively, of shell reference 
surface 

shell stress resultants 

ring hoop force 

ring potential energy 

shell displacements in meridional, circumferential, 
and normal directions, respectively 

ring centroidal displacements 

ring rotations 

effective shell force loads per unit of area 

axial and radial coordinates , respectively 

normal distance of reference surface from shell 
inner surface 

shell stretching strains 

ring hoop strain 
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n (i) n (i) n (D 
01 ,02 ,©12 

9 1 » & 2> e 12 
e 


St 


e <t> 

K ,K ,T 

x’ y 

K 1 > K 2 > K 1 2 
X 

*0 

X* 




i] 

Vi ,v 2 


5»n,v 


p 

P S » a (()» a Z 
°S<j) >°sz > a <pz 

x.'I'.e 

0) 

Vectors: 

F 


effective thermal loads (i = 0 or 1), eqs. (9) 

effective free thermal strains 

stringer free thermal strain 

ring effective free thermal strain 

ring bending strains 

shell bending strains 

load factor (for proportional loading) 

load factor for nonlinear prebuckling state 

limit load 

orthotropic shell wall normal stiffness coefficients 

X - X 0 

eigenvalues 

orthotropic shell wall shear stiffness coefficients 

Poisson contraction ratios with meridional or 
circumferential stress acting, respectively 

shell displacements in axial, radial, and 
circumferential directions, respectively 

mass density 

three-dimensional normal stress components 
three-dimensional shear stress components 

shell rotations about circumferential, meridional, 
and normal directions, respectively 

frequency of harmonic vibrations 
Y' 
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i 


eight (or six) element column vector of dependent 
variables {P,Q,S,Mi ,£,n,v,x> 

particular solution 

complementary solutions 

4x4 (or 3x3) Matrices: 

[B] j [D] boundary condition matrices, eq. (25) 

[D] effective [D] for interior boundaries, eq. (48a) 

[D] additional effective [D] for closed branch 

boundaries, eq. (55) 

ring eccentricity matrix, eq. (34a) 

ring stiffness matrix, eq. (30a) 

matrices relating {c}, { d } of first subinterval 
to that of final sub interval of a closed branch, 
eqs. (59) 

[S] scaling matrix for supplemental conditions, 

eqs. (42) 

[U] ,[W] force and displacement submatrices of 

complementary solutions, eq. (36a) 

[V] ,[Z] additional force and displacement submatrices of 

complementary solutions required on closed branch 

[k] ring prestress matrix, eq. (103) 

[y] ring mass matrix, eq. (112) 

4x1 (or 3x1) Matrices: 

{ c } superposition constants, eqs. (37) 

{d} additional superposition constants for a closed 

branch, eqs. (52) 

{G},{J} force and displacement submatrices of particular 

solution vector, eq. (36b) 


[e] 

[k] 

[p]»[p] 
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{L} 

{£> 

is,) 

e 

U f } 


effective boundary loads, eq. (25) 

effective (L) for interior boundaries, eq. (48b) 

nonhomogeneous ring matrix due to ring 
eccentricity and thermal loads, eq. (34b) 

nonhomogeneous ring matrix due to mechanical 
loads, eq. (30c) 

nonhomogeneous ring matrix due to thermal loads, 
eq. (30d) 


<i) 


{& associated with externally applied loads 

nonhomogeneous matrix relating {c} of first sub- 
interval to that of final subinterval of a 
closed branch, eq. (59a) 


ring displacements (u^u^.u^ ,w^} 


{u} 

{y} shell forces {P,Q,S,Mi} 

{z} shell displacements (5»n,v,x) 

Generalized field variables and operators: 


H(e) 

Lj (u) 

L 2 (u) 

L ll(u»v) 

9l(u) 

u 

e 

a 


V e i’ a i 


linear operator relating stresses and strains 

linear operator representing linear part of the 
strain-displacement relations 

quadratic operator representing the nonlinear part 
of the strain-displacement relations 

bilinear operator defined by the identity 
L 2 (u + v) “ L 2 (u) + 2Ln(u,v) + L 2 (v) 

linear operator representing conservative live loads 

displacement 

strain 

stress 

eigenfunctions 
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Subscripts : 


0 

prebuckling state variable 

1,2,3 

meridional, circumferential, and normal components, 
respectively (same as s,<j>,z) 

<V> 

estimate after k iterations 

Superscripts : 


( > (a) 

antisymmetric component 

< > (k) 

3 k ( )/3X k 

( ) <s) 

symmetric component 

( > T 

transpose 

(") 

load or linear response variable at unit X 

( )' 

3( )/3s 

• 

( ) 

3( )/3<t> 

• 

( ) 

3 ( )/3r 

Matrix subscripts : 


0 

evaluated at the initial point of a subinterval 

1 

evaluated at the final point of a sub interval 


GOVERNING EQUATIONS 


Mathematically speaking, elastic response problems of shell structures 
are boundary-value problems in differential equations. In general, to 
formulate such problems, it is necessary to start with a geometrically 
nonlinear shell theory, i.e., one valid for rotations of moderate size.* 
Also, an analogous theory for elastic rings must be available to formulate 
boundary conditions associated with ring attachments. As a preliminary to 
the formulation of specific types of response problems solved by the SRA 
programs, suitable nonlinear theories of shells of revolution and rings 
are presented in this section. 


*In this approximation, both the strains and rotations are small compared 
to unity, but the rotations may considerably exceed the strains. 


8 


Shell Equations 


Nonlinear shell theories have been developed by Sanders (ref. 14) 
and others. However, for the purpose of numerical analysis of shells 
of revolution, it has been shown that Novozhilov's shell equations (ref. 

15) have the advantage that, by the proper choice of dependent variables, 
explicit reference to the meridional radius of curvature can be elim- 
inated (ref. 9). 

In reference 11, Novozhilov's equations have been generalized, 
through the principal of virtual work, to include the nonlinear case of 
moderate rotations. For numerical analysis, it is convenient to trans- 
form the equilibrium and kinematic equations into a set of eight dif- 
ferential equations in eight basic force and displacement shell variables 
referred to fixed coordinate directions. Four of these variables are the 
effective shell forces in axial, radial, and circumferential directions, 
denoted as P, Q, and S respectively, and the meridional bending moment 
Mi. These components act on normal sections tangent to small circles of 
the shell reference surface (fig. 1) and are all measured per unit of 
circumferential length along the small circle. The remaining four 
variables are the analogous reference surface displacements, denoted as 
£, n, and v, and rotation x» These variables, as well as the notation 
used for other shell variables, are shown in figure 2. As shown, s,<f> 
reference surface coordinates are used where s measures meridional arc 
distance from a reference small circle and <t> measures circumferential 
angle from a reference meridian. The normal distance z measured from the 
reference surface completes the three-dimensional triad of directions. 

The transformation of the equations is accomplished with the use of 
the Gauss-Codazzi surface compatibility relations. Employing the prime 
and dot to denote partial derivatives with respect to s and <j>, respectively, 
the resulting equilibrium equations are 

(rP)' + (r/R 2 )S* - (2/r)M 12 * - (r'/r)M 2 "' - (r/R 2 ) [ (T 3 + T 2 )6]‘ 

+ r'(T 2 t|» + T 12 x) * + r[(r/R 2 )X! - r'X 3 ] + r'Li* = 0 

(rQ)' + r'S* - T 2 + M 2 "/R 2 - r'[(Ti + T 2 )6]’ - (r/R 2 )(T 2 i|> + T 12 x)* 

+ r[r'Xi + (r/R 2 )X 3 ] - (r/R 2 )L!* = 0 (1) 

(rS)' + r'S + T 2 * + M 2 */R 2 - r'(T! + T 2 )e - (r/R 2 ) (T 2 ^ + T 12X ) 

+ rX 2 - (r/R 2 )L! = 0 

(rMi)' + r[r'P - (r/R 2 )Q] - r'M 2 + 2 Mi 2 * - r(T lX + T 12 i|/) + rL 2 = 0 

where the surface force (X 3 ,X 2 ,X 3 ) and moment (L 3 ,L 2 ) components are 
referred to undeformed coordinate directions (fig. 2). 
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The nonlinear terms in equations (1) can be conveniently thought of 
as the following additional load terms applied to the linearized equations. 


Xi = -[ (Tj + T 2 ) 6] ' /r 

(2a) 

X 2 - -(r'/r) (Ti + T 2 )6 

(2b) 

x 3 - 0 

(2c) 

Lj ■ T 2 iJi + T 32 x 

(2d) 

L 2 - -(T lX + T 12 *) 

(2e) 


Additional effective surface loads dependent on the shell deformation 
arise in the cases of an elastic foundation attached to the shell wall and 
loading by a live normal pressure field. An orthotropic elastic 
foundation is considered under the assumption that it produces reactions 
per unit of surface area in meridional, circumferential, and normal 
directions which are proportional to the corresponding shell displacements 
at the surface to which it is attached. It is assumed that the 
attachment surface is the shell inner surface (i.e., the surface of 
inward pointing positive z”direction) . In terms of the displacement of the 
reference surface the foundation loads are 


Xi - 

-ki(u - 

■ *X) 

(3a) 

X 2 “ 

-k 2 (v - 

1 

N > 
-6- 
N— / 

(3b) 

X 3 - 

-k 3 w 


(3c) 

hi = 

zX 2 


(3d) 

l 2 ■ 

-zXi 


(3e) 


where kj, k2, and k3 are foundation moduli and z is the normal distance of 
the reference surface from the inner surface. The effective loads of a live 
pressure field Ap(s,<f>,z), assumed to act at the reference surface, are 


X! - Ap X (4a) 
X 2 - \p\p (4b) 
X3 * A[p(ej + e2) + u3p/3s + v3p/3<|> + w3p/3z) (4c) 
Li - 0 (4d) 
L 2 - 0 (4e) 


where ej and e 2 are the linearized stretching strains in meridional and 
circumferential directions, respectively. 
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Equations (2), (3), and (4) isolate all terms of equations (1) other than 

standard terms of a linear shell statics problem. 

The four basic kinematic equations may be written in the form 


V = r'x + (r/R 2 )ei 

n' = -(r/R 2 )x + r'ei 

v' - -£*/R 2 - (r'/r)(n* - v) + e 12 


(5) 


x' ■ < 1 

where e\ 2 is the linearized shearing strain and <i is the meridional 
bending strain. 

Equations (1) and (5) are eight partial differential equations in the 
eight response variables P,Q,S,Mi ,£,n,v,x. Supplemental equations are 
necessary to express the excess variables of these equations in terms of 
the eight basic variables. The nonlinear strain- rotation equations and 
the partially inverted constitutive equations provide some of the 
supplemental equations. These are 


ei - ei - (1/2) (x 2 + 6 2 ) 

e 2 - e 2 - (1/2) (if; 2 + 0 2 ) (6) 

e 12 = e 12 " X'P 

and 

T 2 = XhE 2 + ^ 12 k 2 + *13^1 + *14^1 “ © 2 ^ (7a) 

M 2 = A 2 2£ 2 ■*" ^22 K 2 ■*" ^23^1 ^24^1 ~ ®2^^ (7b) 

e l = ^31 e 2 + ^32 k 2 + ^ 33^1 + ^34 M 1 (7c) 

— X44G2 "t" ^42 k 2 ■*“ ^43^-1 + ^-44^1 (7d) 

2Mi 2 = y 1 1 k 2 2 + yi 2 S - 20 i 2 ^ 1 ^ (7e) 

G 12 = y21 K 12 + V*22S (7f) 


where ei,e 2 ,ei 2 and ki,k 2 ,k 12 are reference surface stretching and 
bending strains, respectively, and 
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Tl - Ti + 0i 


Mi - Mi + 0i 


( 0 ) 

(1) 


<12 = <1 2 “ g 12/ r 2 

S = S - (1/2) (Ti + T 2 )0 + 0i2 (O) + 20 12 (1) /R2 


(8a) 

(8b) 

(8c) 

(8d) 


The thermal loads 0]/ m ), 0 2 ^ m \ and ©i2^ m \ for n ■ 0 or 1, are given in 
terms of the free thermal strains 0j, 0 2 , and 0j 2 by* 


© 1 ^ = / tE i / (1 - ViV 2 )](0i + v 2 0 2 )z m dz + 0 gt / m ^ 

02^ m ^ = /(E 2 /(l - viv 2 )](0 2 + vi0i)z m dz (9) 

©12 (m) = /E 12 0i 2 z m dz 


where Ej, E2 , Ei 2 , vi , and v 2 are orthotropic shell wall elastic moduli, 
the integrals are through the shell wall thickness, and 0 gt ( m ) are given 
in terms of the stringer free thermal strain 0 st by 


0 „ (0) = (NEA/ 2irr) 0 
st st 


0 (1) = e 0 (0) 

st z st 


(9a) 


Here NEA is the total stringer stretching stiffness and e z is the normal 
eccentricity of stringer section centroids relative to the shell refer- 
ence surface. The stiffness coefficients A^j and in the constitutive 
equations (7) are defined in Appendix A. J 


In addition to equations (6) and (7), the following equations 
complete the supplemental equations. 

T 12 = S - 2M 12 /R 2 - (1/2) (Tj + T 2 )e 
and T 2 = (r/R 2 )P + r’Q 

e 2 * (r) + v*)/r 

K 2 = [r'(x + 5 *'/r) + (v* - n '‘)/R 2 ]/r 
<12 - ®i 2 /R 2 " <X ~ 5/r)*/r 

ip ■ (v - n*)/R 2 + (r'/r)^* 

0 - (r'/r) (v - n‘) - £*/R 2 
u - (r/R 2 K + r’n 
w = -r'5 + (r/R 2 )n 


( 10 ) 


(ID 


* Although the shearing free thermal strain 0 12 is zero for an orthotropic 
material, it will be seen to be convenient to include it in the formulation. 


12 



It may be noted that the only nonlinear terms appearing in equations 
(5) through (11) are those in equations (6),(8d), and (10). Just as the 
nonlinear terms in the equilibrium equations (1) may be viewed as 
additional mechanical loads applied to the linearized equations, the 
nonlinear terms in equations (6) and (8d) may be viewed as additional 
thermal loads applied to linearized versions of these equations. Noting 
that the three-dimensional strains q, £2> e 12 appear in the stress-strain 
relations as ei - 0i, e 2 ~ 0i>and tj2 ~ ®12> it follows that the nonlinear 
terms of equations (6) are equivalent to the following additional free 
thermal strains applied to the linearized equations, 

01 = -d/2) (x 2 + e 2 ) 

e 2 = -d/2) (^ 2 + e 2 ) (12a) 

012 - -x* 

which do not vary through the shell thickness. The corresponding thermal 
loads are obtained by substituting equations (12a) into equations (9) to give 

©1 (m) = -(l/ 2 )[C s (m)( x 2 + 0 2) + c 12 ( m )d 2 + 02)] 

G 2 (m) . _( 1 / 2 )[c 2 (m)(^2 + e 2) + c 12 ( m )(x 2 + e 2 )] (12b) 

01 2 (m) = -G 8 (“> X i|» 

where the stiffness coefficients Cg( m ), C 12 ( m ), c 2^ m ^» and G g ( m ) are 
defined by equations (A-3) of Appendix A. The nonlinear term in equation 
(8d) is evidently equivalent to the additional thermal load 

012 (O) = -(1/2) (Tj + T 2 )0 (13) 

Finally, to identify the nonlinear term in equation (10) as an 
additional load, note that the only places where T 12 is required are in 
equations (2d) and (2e). It therefore follows that this nonlinear term is 
equivalent to the additional loads* 


Li = -(1/2) (T X + T 2 )x6 

1*2 - (1/2) (Ti + T 2 )*6 


(14) 


Thus, the general nonlinear field equations may be viewed as a standard 
set of linearized equations with the additional load terms given by 
equations (2), (3), (4), (12b), (13), and (14). 


*Sinca these loads are smaller then similar terms in equations (2d) and 
(2e) by a factor of the rotation 6 , it is consistent with the moderate 
rotation theory to neglect the nonlinearity in equation (10). 
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Ring Equations 


When ring stiffeners are attached to the shell, boundary conditions 
for the shell equations must be generated which represent the ring 
behavior. Ideally, the ring reactions enter the shell at a single 
meridional station, the ring boundary, on the shell reference surface 
at which the shell displacements are continuous and related to the shell 
force jumps in accordance with the governing ring equations. A set of 
suitable ring equations are derived in this section. 

As are the shell equations of the previous section, the ring 
equations are based on moderate rotations and are derived through a 
principal of virtual work. These equations are based on the following 
assumptions . 

(1) All geometrical and mechanical properties of the ring are 
axisymmetric. 

(2) The ring material is homogeneous and isotropic. 

(3) The effects of nonuniform warping of ring sections, trans- 
verse shear strains, and shear center eccentricity relative 
to the section centroid are neglected. 

The origin of ring cross-sectional x,y-axes is supposed to be at the 
centroid of the ring section, i.e., /xdA ■ /ydA * 0 where A is the section 
area. With respect to right-handed shell coordinates s,<|>,z, x is chosen 
positive in the axial direction acute to the positive (or negative) s- 
direction if the positive z-direction points away from (or towards) the 
axis of revolution, and y is chosen positive in the radial direction 
pointing away from the axis of revolution (see fig. 3). 

For a one-dimensional theory of rings the centroidal hoop strain 
e<j, is the only stretching strain of consequence. The strain-displacement 
relations are 


e <|> = + u y )/ a + (1/2) (w x 2 + w y 2) 

K x " ( u y ~ a ^)/a 2 
Ky = -(u"/a + wp/a 
x - (w^ - u x /a)/a 


(15) 


where k x and Ky are the bending strains of the centroidal axis, in and 
out of the plane of this axis respectively, and x is the twist per unit 
of circumferential length. Neglecting transverse shear strains, the 
rotations w x ,w y may be written in terms of displacements as 
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(u y - u ( j ) )/a 


( 16 ) 


w y - -u^/a 

Integrating by parts the following expression for the virtual change in 
potential energy 

.2 ir 

6U - So ( T <J> 6e <j, + Mx«*x + My6iCy + M^x)ad<(. 

,2ir 

"Jo < F x 6u x + V u y + + N x 6w x + N y 6w y + N^w^ad* (17) 

and applying the principle of virtual work, 6U = 0, yields the following 
equations expressing equilibrium of forces and moments in the undeformed 
coordinate directions (fig. 3). 

(M^/a - M^/a + N y - w y T^)' + aF x = 0 
- N x + w x T\)* - Tj, + aF y ■ 0 

(18) 

tJ - M^/a - N x + w x T + + aF<, - 0 
Mj + My + aN^ = 0 

For a one-dimensional ring theory, the constitutive equations are 
unchanged from those for a straight elastic bar. Neglecting the effect 
of nonuniform torsion, for a homogeneous, isotropic bar these are (ref. 16) 


t<j, - EA( e<j) - ep 

Mx = EI x K x - EI X yKy 
My S — EI x yK x t ElyKy 

M^ = GJt 


(19) 


where the ring free thermal strain 0^ is assumed to be uniform over each 
cross section. 

In analogy with the nonlinear shell equations, the nonlinear terms 
in equations (15) and (18) may be viewed as effective additional mo- 
ments and free thermal strain applied to the linearized ring equations. 
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These additional loads are 


N * 
x 


-W T 
X <p 


N - - w T . 

y y 4> 

EAd^ - -(1/2) EA(w x 2 + w y 2 ) 


(20a) 

(20b) 

(20c) 


To reduce the ring equations to a more useful form, equations (15) , 
(16) , and (19) are substituted into equations (18) to eliminate all 
response variables except u x , Uy, , and . With the understanding that 
the nonlinear terms are represented in the load terms according to 
equations (20) , the resulting equations are 


• • •• • 

(El u“ - GJu ) + El (u“ - u") + a (El + GJ)w" = a 3 (aF + N*) 

y x x xy y <t> y 4> x y 

El (u" + aw") + a 2 EAu + El u" + a 2 EAu A - El u" 

xy x <j> y xy <J> x <t> 

• a 3 (aF - N* + EA6 a ) 

y x <f>' 

EI xy (u x + aw ^ + (EI x U y ' a2EAu y> - ( a2EA + EI X > U *‘ 

= a3 < aF 4, - N x - EA %) 

(Ely + ^“x + EI xy (u y " V + a(EI y W <j> ~ GJ V = a3N <^ 


SOLUTION OF EQUATIONS 


In this section the governing equations of the previous section are 
specialized for the different modes of response treated, and the corres- 
ponding methods of solution are presented. 


Linear Asymmetric Response (SRA 100) 


This program solves linearized versions of the shell and ring equations 

subject to harmonic mechanical and thermal loads. Since all load and 

response variables are periodic functions of a 4> with period 2tt, they may 

be represented as Fourier series in the form ?_(A cos nd> + B sin ni>), 

n*U n n 

where the harmonic amplitudes A and B are in general functions of s. 

n n 
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If the loads have an axial plane of symmetry, say <P = 0, then the 
Fourier series for each load term reduces to a sine or cosine series. 

A symmetrical loading is defined as one for which the expansions for Xi , 
x 3» l 2»01 » e 2»F jN , ,0, (denoted henceforth as normal type load variables) 
are cosine seriel, whereas the expansions for the remaining loads X2,Li,8i2 
F^.Nx.Ny (denoted henceforth as shear type load variables) are sine series.* 
The reverse is true in the case of an antisymmetrical loading, and a gen- 
eral load consists of both symmetric and antisymmetric components (table I). 

A symmetrical response is defined as one for which the expansions for 
P,Q,Mi,£>ri>X>vix>Uy> w ({> (denoted henceforth as normal type response variables) 
are cosine series, whereas the expansions for S,v, and u^ (denoted hence- 
forth as shear type response variables) are sine series. The reverse is 
true in the case of an antisymmetric response, and a general response 
consists of both symmetric and antisymmetric components (table I). 

Inspection of the nonlinear shell and ring equations shows that a 
symmetric loading gives rise only to a symmetric response. For the 
linearized equations, it is also true that an antisymmetric loading gives 
rise only to an antisymmetric response. Furthermore, for the linearized 
equations , the response to each load harmonic is a pure harmonic of the 
same wave number. 

Symmetric load- response equations .- In this section the boundary- 
value problem for the symmetrical response to the n-th harmonic of the 
symmetrical load components is formulated. For the sake of simplicity, 
the same symbols as used previously for physical load and response 
variables will be used to denote the corresponding harmonic amplitudes. 

Differential equations: Substitution of the symmetric load and 

response components for the n-th harmonic into the linearized form of 
equations (1) and (5) gives the following ordinary differential 
equations. 


(rP)' + n(r/R 2 )S + n 2 (r'/r)M 2 - 2(n/r)Mi2 + r[(r/R 2 )Xi - r'Xa] 

+ nr'Li * 0 

(rQ)' + nr’S - (n 2 /R 2 )M 2 - T 2 + r[r'X! + (r/R 2 )X 3 ] - n(r/R 2 )L! =0 

(22a) 

(rS)’ + r'S - (n/R 2 )M 2 - nT 2 + rX 2 - (r/R 2 )L! = 0 
(rMj) ' + r[r'P - (r/R 2 )Q] - r'M 2 + 2nMi 2 + rL 2 - 0 


*It will be convenient in the remaining discussion to refer to free thermal 
strains 0i,02»0i2> and e <j> as loads. More precisely, the thermal loads are 
given in terms of the free themal strains by equations (9) for the shell 
and as EAO^ for rings. 
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V - r'x - (r/R 2 )ei = 0 
n' + (r/R 2 )x - r'ei = 0 
v’ - (n/R 2 )C - (r'/r)(v + nn) - ei 2 = 0 

X r - =0 


(22b) 


In these linearized equations, the elastic foundation loads, equations 
(3)» are considered to be included in the load terms with the Riven 
applies loads; however, the live load terms, equations (4), are 
neglected. 

The supplemental equations, expressing the excess variables of 
equations (22) and (3) in terms of the eight basic variables, are 
equations (6), (7), and (8) with nonlinear terms omitted, plus the 
following from equations (11) 


u = (r/R 2 )£ + r'n 

(23a) 

w - -r'£ + (r/R 2 )n 

(23b) 

Ti - (r/R 2 )P + r'Q 

(23c) 

re 2 = n + nv 

(23d) 

ri< 2 = r'(x - n 2 £/r) + n(nn + v)/R 2 

(23e) 

r(Ki 2 - e 12 /R 2 ) = n(£/r ~x) 

(23f) 

\p = (nn + v)/R 2 - nr'5/r 

(23g) 


Equations (22) and the supplemental equations are a system of eight 
first-order differential equations which may be written compactly in 
vector form as 


Y' - F(s,Y) (24) 

where Y is the eight-element column vector (P,Q,S,Mi ,£ ,n ,v,x) [fig. 2(a)]. 

Boundary conditions: Branch edges, branch points, the closure point 

of a closed branch, and the location of interior rings or other meridional 
discontinuities are defined as boundaries. Additional artificial points 
of subdivision of the meridian may be required to limit subinterval length 
so that the small difference of large numbers does not occur in the 
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superposition of complementary and particular solutions of equations (24) 
(see p. 22). In general, the region of integration of equations (24) 
consists of a main branch and subsidiary branches. The main branch is a 
continuous line consisting of segments of the shell reference meridian 
which in the case of only open branches begins at some arbitrary edge 
a xid terminates at some other arbitrary edge. If the meridian contains a 
closed branch (only one is allowed) , the closed branch is the main 
branch, which begins at some arbitrary nonbranching point and terminates 
at the same point. At a branch point, only one branch exits the branch 
point, i.e., has increasing s-values away from the branch point. All 
other branches intersecting the branch point must enter the branch point, 
i.e., have s-values increasing towards the branch point. All branches 
entering a branch point are described by increasing s before the exiting 
branch. 

General linear boundary conditions for each boundary may be written 
in the form 

[B J A { y } + [D] {z} = {L> (25) 

where (y } and {z} are 4x1 force and displacement subvectors of Y and 

A{y} = ±{y} at edges (26a) 

A{y} = (y} + - £ Cy) at interior boundaries (26b) 


In equation (26a) , the minus sign applies only at the terminal edge (if 
one exists) of the. main branch; in equation (26b), (y} + is the value of 
(y } at the boundary on the exiting branch, and E{y} - is the sum of the 
values of y on the branches entering the boundary. As implied by the form 
of equation (25), at interior boundaries the displacement vector {z} is 
continuous, i.e. 

(z} + = {z}' (27) 

The matrices [B], [D], and {L} are generated by SRA 100 in the case of 
force-free or ring boundaries or dome closure edges. The first two types 
are discussed here, whereas dome closures are discussed in Appendix B. 

At force-free boundaries A{y} = {0}, so that 

[Bj = [I] 

[D] = [0] (28) 

(L> = {0} 
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Substitution of the symmetric load and response components for the 
n-th harmonic into the linearized ring equations (21) yields 


[k]{u> - U f } + U t } 


( 29 ) 


where [k] is the ring stiffness matrix given by 


[k] 


n 2 (n 2 EI y +GJ)/a 2 


1 

a 


Symmetric 


n 4 EI /a 2 n 3 EI /a 2 

xy xy 

EA+n 4 EI /a 2 n(EA+n 2 EI /a 2 ) 

X X 

n 2 (EA+EI x /a 2 ) 


-n 2 (EI y +GJ)/a 

-n 2 EI /a 
xy 


-nEI / i 
xy 


El + n 2 GJ 

y 


(30a) 


and (u) , {if} and (£ t ) are ring centroidal displacement, and mechanical 
and thermal load vectors given by (see fig. 3) * 


{u} 


u 


< 7 > 

u +f 

l w W 


u £ } 



U t ) 


EA6. < 

<P 


0 

1 

n 

V°V 


(30b) 


(30c) 


(30d) 
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In order to derive boundary conditions of the form of equation (25) 
from equation (29) , it is necessary to relate the ring centroidal dis- 
placement and load vectors, {u} and {if}, to the shell reference surface 
displacement and force jump vectors, {z} and My}, at the corresponding 
boundary point. Equilibrium of forces and moments at the ring centroid 
gives, in terms of the ring eccentricities e x and &y [fig. 3(a)], 


U f } * [B]My} + U f (e) } 


(31) 


where 


a 


[B] 



0 

0 


0 

a 

0 


-ne 0 

x 

-ne 0 

y 

r 0 


(32) 


L -ae ae 0 a 

$ X J 

and {£, } is given by equation (30c) with the ring forces and moments 

per unit of circumferential length replaced by corresponding external 
forces and moments . Furthermore, assuming that the ring centroid is 
connected to the corresponding boundary point on the shell reference 
meridian by a rigid link with the ring free thermal strain 0, , one 
obtains the kinematic relation ^ 


{u} - [e] {z} + U } 
e 

where [e] is an eccentricity matrix given by 


(33) 


[e] 


1 

0 

ne /r 
x 


0 

1 

ne /r 

y 

0 


0 

0 

a/r 

0 


-e 

a 

0 

1 


{ V 


0 , < 


(34a) 


(34b) 
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Substitution of equations (31) and (33) into equation (29) yields 
the desired boundary condition in the form of equation (25), where [B] is 
given by equation (32) and 

[D]--[k][e] (35a) 

{L} =» -U, (e) } - U. } + [k]U } (35b) 

r t e 

As a check of this result, it is noted that [B]”^ = [e]^/r. Multiplication 
of equation (25) by [Bj“l shows that the self-adjointness (i.e. symmetry) 
of equations (29) is preserved. 

Method of solution .- The original method of solution of the linear 
boundary-value problem, equations (24), (25) and (27), is denoted here as 
the Gaussian elimination method (ref. 9). This method consists of 
subdividing the range of integration (i.e., the shell meridian) into a 
number of suitably small subintervals, the end points of which have been 
denoted as boundaries in the previous section. A forward integration 
scheme, such as Runge-Kutta, is used to integrate equation (24) over each 
subinterval between consecutive boundaries, to obtain eight linearly 
independent complementary solution vectors Y c ^, k = 1,...,8, and a 
particular solution vector Y . Initially, at the starting point of each 
subinterval, the matrix of column vectors (Y c ^)] is chosen to be the 8x8 
identity matrix and Yp the 8x1 null matrix. The boundary conditions (25) 
and (27) are used to set up a system of algebraic equations for the 
constants of superposition for each subinterval. These are solved 
efficiently by Gaussian elimination in terms of 4x4 matrices, and the 
results used to superpose the complementary and particular solutions to 
obtain the desired solution. The sub intervals must be small enough so 
that the superposition of solutions <bes not involve taking the small 
differences of large numbers, with a consequent loss of significance. It 
is characteristic of this method that the information contained in the 
boundary conditions is not used during the forward integration of the 
differential equations and that the initial conditions used for the 
complementary and particular solutions are arbitrary within the condition 
that the 8x8 initial value matrix [Y^k)] should be nonsingular. 

Later, Zarghamee and Robinson (ref. 7 ) proposed the use of starting 
conditions for the complementary and particular solutions which imply 
satisfaction of the boundary conditions. Since four conditions are 
known at the initial edge, they reasoned that only four complementary 
solutions are required to satisfy the four conditions at the final edge. 

On the other hand, since only four conditions are known in terms of eight 
variables at the initial edge, there is still some arbitrariness in the 
determination of the starting conditions in this method.* Their technique 


*As a consequence, the problem of "long sub intervals" noted above remains 
in the Zarghamee method. It is noted here that a new method, termed the 
field method (ref. 17), which eliminates this problem as well as providing 
other benefits, is currently being investigated. 
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was generalized to general linear boundary conditions and open branched 
shells by Anderson, et al. (ref. 2, Appendix A). During the course of the 
present study, it was found that the supplemental starting conditions 
proposed in reference 2 lead to the inversion of poorly conditioned 
matrices. In this section the Zarghamee method is presented with new 
supplemental starting conditions, and it is also generalized to include 
shells with a closed branch. 

Open branches: The 8x4 matrix of the four complementary solution 

vectors [Y C W], k - 1,...,4, and the 8x1 particular solution matrix Yp 
are partitioned into 4x4 submatrices U,W and. 4x1 submatrices G,J as 
shown below (for simplicity in writing, the brackets and braces used for 
4x4 and 4x1 matrices, respectively, are omitted in the remainder of this 
section) . 


t?c (k)) ■ [- 5 -] (36a) 

Y P - {-f } (36b) 

Then the desired solution y,z may be written for the i-th subinterval as 


y * G + Uc ± (37a) 

z - J + Wc i (37b) 

where c^ is a 4x1 matrix of superposition constants for the i-th 
sub interval. 

At a starting edge of an open branch, the boundary conditions (25) 
may be written 


Byo + Dzq = L (38) 

where the subscript 0 denotes initial values. Substitution of equations 
(37) into equation (38) shows that equation (38) will be satisfied 
regardless of the value of c for the subinterval considered if 

BU 0 + DW 0 * 0 (39a) 

BG 0 + DJ 0 = L (39b) 
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Equations (39) are then starting conditions for the matrices, U,W,G,J, at 
an initial edge. 

Equations (39a) and (39b) are respectively 16 equations in 32 
unknowns and 4 equations in 8 unknowns and hence do not have a unique 
solution. In order to formalize the procedure, it is necessary to augment 
equations (39) with supplementary conditions such that 

(1) the initial values Uo,W 0 ,G 0 ,Jo are uniquely determined, and 

(k) 

(2) the complementary solution vectors are linearly independent. 

In Appendix C, it is shown that condition 2 will be satisfied for any 
supplementary condition for Uo,Wq of the form ciUq + $Wq * I, where a and 
6 are 4x4 matrices. In order to minimize the calculation of the initial 
values, for any particular choice of a and 8, it will be convenient to 
choose the supplemental conditions for Go, Jo as aGo + 8Jo * 0. If B is 
nonsingular, as is the case if no kinematic constraints are specified at 
the boundary, a suitable set of supplemental conditions are obtained by 
simply taking a ■ 0 and 8*1, viz. 


W 0 - I (40a) 

Jo - 0 (40b) 

Substitution of equations (40) into equations (39) gives 

U 0 * -B -1 D (41a) 

G 0 * B -1 L (41b) 


If B is singular, equations (40) are replaced by 

±U 0 + SW 0 - I (42a) 

±G 0 + SJ 0 - 0 (42b) 

where S is a diagonal scaling matrix, the purpose of which is to provide 
dimensional homogeniety to equations (42) . The first three diagonal 
elements of S are taken to be Cj'^Vt, and the fourth diagonal element 
is where t is an effective thickness given by 

t - (12C 1 (2) /C 2 (0) ] 1/2 (43) 

and Cj (0) and C*^) are meridional stretching and bending stiffnesses 
[see eq. (A-3) of Appendix A], 
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Substitution of equations (42) into equations (39) gives the initial values 


W 0 = (BS + D) -1 B 

(44a) 

U 0 = ±(I - SW 0 ) 

(44b) 

J 0 = +(BS + D) _1 L 

(44c) 

Gq = +SJ o 

(44d) 


Since, conceivably, either BS + D or BS - D may be singular, both upper 
and lower signs in equations (42) and (44) are allowed. Equations (44) 
could of course be used as initial edge starting values in all cases; 
however, the relative simplicity of equations (40) and (41) suggests their 
use in the common case of nonsingular B. 

At an interior boundary, at which several open branches may intersect, 
the boundary conditions (25) and (27) may be written as (see fig. 4) 

J 

B[y ij+1,0 " ^ y i j ,l 1 + Dz ij+1»0 = L (45a) 


Z ij+1,0 " Z ij,l 


j = 1,***,J (45b) 


where ij are the subinterval numbers (generally nonconsecutive) of sub- 


intervals terminating at the boundary, the number of which is denoted as 
J. In equations (45) the first subscript refers to the subinterval 
number and the second subscript 0 or 1 indicates evaluation at either 
the beginning or end of the subinterval, respectively. Substitution of 
equation (37b) into (45b) gives the c^ for entering subintervals in 

j * * 


terms of Cj_ for the exiting sub interval, viz. 
1 



W 


-1 


1 d » 1 


(J 


ij+1,0 


" J i x + 

V 1 


^j+l.O^j+l^ 


j * 1, • • • , J (46) 


Substitution of equations (37) into equation (45a) and elimination of ^ 
through use of equations (46) , shows that equation (45a) will be J 
satisfied regardless of the value of ^ if 


*It is shown in Appendix C that W - are nonsingular matrices. 

1 j * 1 
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V 1 ’ 0 


(47a) 


where 


BU i J+ l,0 + DW i J+ l,0 " ° 

BG i J +l,0 + D J i j+1 , 0 ■ L 


D « D - b£ A. 

J-1 i 

J 

L = L + B^ A. 

j“l J 


and 


\ ' VV 


-1 


A. " G. I “ A. J. , 

*j V i j» 1 


(47b) 

(48a) 

(48b) 

(49a) 

(49b) 


Equations (47) are the starting conditions for the matrices U,W,G,J on the 
exiting branch of an interior boundary. Since they are of the same form 
as equations (39) , the starting values of these matrices are also given by 
equations (40) and (41), or (44), with D and L replaced by 6 and L, 
respectively.* 


If the shell contains no closedbranch, a terminal edge will be 
reached at the end of, say, the mrth sub interval. For this boundary, the 
boundary condition (25) may be written as 



+ Dz 


m,l 


L 


(50) 


*It may be noted here that since the starting conditions for the 
complementary solution matrices U,W are independent of the boundary load 
vectors L, U and W are independent of all load (nonhomogeneous) terms in 
both differential equations and boundary conditions. (The same will be 
seen to be true for the additional complementary solution matrices V,Z 
required on closed branches). Consequently, in a sequence of problems 
in which only load terms change, the complementary solutions need be 
computed just once. 
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Substituting equations (37) into equation (50) and solving for gives 


c 

m 


(DW m , - BU ,)" 1 (L - DJ , + BG ,) 
m, x m j X- m y x m y x 


(51) 


Starting with this value of c m , equation (46) is used recursively to 
obtain the for each subinterval, after which the solution is given 
by equations (37) . 

Closed branches: The present method requires the calculation of four 

additional complementary solution vectors on a closed branch. If the 
matrix of these vectors is partitioned into 4x4 submatrices V and Z, the 
desired solution may be written for the i-th subinterval on a closed 
branch as [cf. eqs. (37)] 


y - G + Uc ± + Vd ± (52a) 

2 = J + Wc ± + Zd i (52b) 


where d^ is an additional 4x1 matrix of superposition constants for the 
i-th sub interval. 


Equations (45) are the proper boundary conditions for an interior 
boundary of closed branch. Since only one closed branch is allowed and 
this must be chosen as the main branch (see p. 19), it follows that in 
this case, subintervals ii and ij+1 are entering and exiting closed 
branch subinf ervals , and ij , j = 2,***,J, are entering open branch 
subintervals. Substitution of equations (52) into equations (45) shows 
that equations (45) will be identically satisfied with respect to 

(j = 1,***,J+1), d-j^, and d i +i if c^ is given by equation (46) 

with j ■ 1, equations (47) are satisfied, and in addition* 


d± i = ^r 1 1 Zi j +l ’° dl j +1 

c i. " W i,,l 1(J i +1,0 " J i.,l + ^ij.+l,0 c ij+l + Z i J +l,0 d i J +l ) 


j j 


j - 2, * * * »J 


BV. . . + DZ. .. A = 0 

iyl-1,0 Xj+1,0 


(53a) 

(53b) 

(54) 


*It is shown in Appendix C that Z. . are nonsingular matrices. 

V 1 
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where 


D-D - B(V, n Z “ 1 + l A ) (55) 

V A V 1 j-2 x j 

In Appendix C it is shown that eight linearly independent complementary 
solution vectors (in subinterval ij+1) satisfying equations (47a) and (54) 
do not exist if the corresponding boundary condition matrix B is singular. 
Therefore, singular B matrices (i.e., kinematric constraint) are not 
allowed at boundaries on a closed loop except at the closure (terminal) 
point. Furthermore, in Appendix C it is shown that the linear independence 
of the eight complementary solution vectors Y c (k) depends on the non- 
singularity of the initial values of W and Z. Therefore, the supplemental 
conditions for equations (47) and (54) on a closed branch are always chosen 
to be 


W 

iyKL,0 


Z ij+1,0 


- I 


0 


(56a) 

(56b) 


Substitution of equations (56) into equations (47) and (54) then gives the 
remaining initial values 


U 


ij+i , 0 


V 1 ’ 0 


V 1 ’ 0 


-B~ 1 D 


- 1 " 
-B D 


-B _1 L 


(57) 


The integration on a closed branch is started at an arbitrary 
(nonbranching) point with the initial values 



Z 


1,0 


I 


(58a) 




(58b) 



(58c) 


The solution for these matrices is continued by forward integration from 
the initial point to the final (closure) point. At intervening boundaries 
on the closed branch, the integration is restarted with the initial 
conditions given by equations (56) and (57); at intervening boundaries 
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on open branches, the initial conditions are given by equations (40) and 
(41), or (44) for edge boundaries, and these same equations with D and L 
replaced by D and L for interior boundaries. As each boundary on the 
closed branch is passed, equations (46) with j = 1 and (53a) are applied 
to generate relationships giving cjjd^ in terms of c^d^ of the exiting 
subinterval, viz.* 


c i ' p k c k + \ 


d i " p k d k 


(59a) 

(59b) 


In view of equations (56) , equations (46) for j = 1 and (53a) reduce to 
(setting i^ = k and ij + 1 = k + 1) 


c k = w k,r 1(c k+i • 

d k = Z k,l d k+l 


(60) 


Substitution of equations (60) into equations (59) and comparison with 
equations (59) with k replaced by k + 1 gives the recursion relations 

p k+l = p k W k,l 

A A. __ ^ 

p k+l = P k \,1 (61) 

q k+l = q k “ p k W k,l J k,l 
which are used with the initial values 


Pi " p l * 1 

(62) 

q x = o 

to generate p^, p^, and q^. When the final subinterval is reached, the 

A A 

matrices p = P =. p K’ and q = q R will have been obtained.* 


At the closure point the boundary conditions (25) and (27) may be 


*Here, k is an index for sub intervals on the closed branch only, 
k = 1, • • • ,K. 
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written as Icf. eq. ( 50 )] 


^l.O - + D *i,o ■ L 

z l,0 = z m,l 


( 63 ) 


Substitution of equations (52) and (58) into equations (63) , and 
elimination of and d} through use of equations (59) with k = K gives the 
following two equations for c m and d^. 


B(p - U i ) c + (Dp' - BV )d - L + B(G - q) 
m,l m m ^ 1 in m 5 1 

W -C + (Z . - p)d = -J •, 
m,l m m,l m 

The solution of equations (64) is 

c - (DW , - BU .) -1 (L - DJ + BG ) 
m m,l m,l m,l m,l 

d = (p - Z 1 )' 1 (W m ,c + J ) 

m m,l m,i m m,l 

where 


D = (Dp - BV) (p - Z) 
m,± m,l 


-i 


(64) 


(65a) 

(65b) 





( 66 ) 


It may be noted that equation (65a) for c m is of the same form as that 
for the open branch [eq. (51)]. Starting with these values for Cjn and 
d m , equations (60) are used to obtain and d^ for closed branch sub- 
intervals, and equations (53b) and (46) are used to obtain c^ for open 
branch subintervals. The solution on the closed branch is then given by 
equations (52) and on open branches by equations (37).* 

Antisymmetric loading .- As has been noted on page 17, for the 
linearized shell and ring equations, the response to antisymmetric load 
components is also antisymmetric. It is shown in this section that the 
antisymmetric response can be obtained from the solution of the symmetric 
load-response equations, outlined in the preceding sections. 


*The equations used for the calculation of the components of the three- 
dimensional stress tensor from the solution for the y and z vectors are 
given in Appendix D. 
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Considering as typical normal and shear type loads the normal 
pressure X3 and the circumferential shear X2 , respectively, one has 


X3 = £ IX3 ^ ^ cos n<j> + X3 ^ ^ sin n<|>] 
n=0 

00 

X2 ■ J [X2 ^ sin n<j> + X2 ^ cos n<|>J 
n=0 n n 


( 67 ) 


Here the superscripts (s) and (a) refer to symmetric and antisymmetric 
components, respectively. Considering as typical normal and shear type 
response variables the meridional stress resultant Tt and the shear stress 
resultant T22» respectively, one has 

OO 

Ti - l [ T l n ^ S ^ cos n< f> + T l n ^ sin 

n “° ( 68 ) 

T22 ■ £ [T12 ^ s i n n< f> + T12 ^ cos n$] 

n=0 n n 

From the identities 


sin nd> = cos(n<J> - it/ 2) 
cos n<f> = -sin(n<|> - ir/2) 


(69) 


if follows that the antisymmetric load components are equivalent to 
symmetric components about a rotated plane according to 


X3 ^ sin n<|> = ^ 3 ^^ cos(n<|> - tt/2) 

X2 n ^ 3 ^ cos nd> = ~X2 n ^ a ^ sin(n<(i - ir/2) 

(a) (a) 

Therefore, using the amplitudes X3 *= X3 r and X2 = -x 2 n 

symmetric response equations will give the amplitudes Tj and 
ponding to the solution 

Ti cos(n<t> - tr/2) = Tj sin n<(> 

T12 sin(n4> - tt/2) = -T12 cos n<(> 


(70) 


in the 
T12 corres w 

(71) 
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Comparison of the right hand sides of equations (71) with equations 
(68) shows that 


m ( m 

T l„ - Tl 

T ( a ) _ rp 

T 12 n _t 12 


( 72 ) 


One therefore concludes that the same equations used for the 
symmetric components can be used also to obtain the antisymmetric 
components if the signs of the shear type antisymmetric load amplitudes 
are reversed before the solution and the signs of the shear type anti- 
symmetric response amplitudes are reversed after the solution. Hence, 
the complementary solutions and the [B],[D] boundary matrices for the 
symmetric response components of a particular harmonic can be used also 
for the antisymmetric response components of the same harmonic. 


Buckling of Axisymmetric Equilibrium States (SRA 101) 

This program calculates bifurcation buckling modes of linearized 
asymmetric prebuckling states. The structural loading is assumed to 
have a given spatial distribution, but its magnitude is allowed to vary 
in proportion to a load parameter X. This leads to a set of eigenvalue 
equations for the critical load X Q . A general form of the eigenvalue 
equations for bifurcation buckling and their method of solution are 
presented in Appendix E. In the development presented there, no 
restrictions are placed on the structural geometry or loading, and non- 
linear prebuckling states are included. In this section the iteration 
equations (E-ll) and the inner product [eq. (E-6) , required to calculate 
the eigenvalue estimate according to the Rayleigh quotient, eq. (E-19)] 
are specialized to ring-stiffened shells of revolution assuming linear- 
ized prebuckling states and neglecting prebuckling rotations. 

Iteration equations .- For numerical purposes, the differential form 
of the variational iteration equations (E-ll) is more convenient to use 
and is derived here. First, the differential form of the variational 
eigenvalue equations (E-4) , from which the iteration equations follow, 
is obtained. Applying the usual procedure of taking the difference of 
the governing shell and ring equations evaluated for an initial (pre- 
buckling) equilibrium state and an adjacent (buckling) equilibrium 
state, and linearizing the resulting equations for the perturbation 
variables, leads to the buckling eigenvalue equations. The nonlinear 
shell and ring equations have been presented earlier in the form of 
linear differential equations with nonlinear (and live load) terms 
isolated as additional effective mechanical and thermal loads given by 
equations (2), (4), (12), (13), (14), and (20). Hence, the eigenvalue 
equations are obtained in the same form as the linear system of equations 
[with dead load terms dropped from the linearized form of eqs. (1), (7), (8) 
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and (21)] with effective additional loads derived from equations (2), (4), 
(12), (13), (14), and (20). Since, however, prebuckling rotations are 
neglected in this analysis, the contributions from equations (12), (13), 
(14), and (20c) are not retained. This leads to the following addition- 
al effective loads: from equations (2) 

Xi = -A[(Ti + T 2 )6j7r 
X 2 = -A(r'/r)(T 1 + T 2 )9 

X 3 = 0 (73) 

L 1 = A (T 2 9 + T 12X ) 

L 2 = -A(Tix + T 12 )]j) 


from equations (4) , for live pressure loading. 


Xi = Ap X 
X 2 = Ap9 

X 3 = Ap(ei + e 2 ) 
L i = L 2 = 0 


(74) 


which are identical in form to equations (4) except that pressure gradient 
terms are neglected, and from equations (20a, b) 


N 

x 


-AT.w 

9 x 


N 

y 


-AT.w 
9 y 


(75) 



= 0 


In equations (73), (74), and (75), unbarred variables represent buckling 
mode variables, and barred variables represent unit load prebuckling 
state variables.* 


The iteration equations are then obtained by setting A = 1 in 
equations (73), (74), and (75), interpreting the unbarred variables in 


*Note that since only linearized prebuckling states are treated, in the 
notation of Appendix E, Aq =0, uq = oq = 0, and A = y. 
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these equations as being known inputs from the previous [(k - l)-th] 
iteration, and solving the linear system of equations with these loads for 
the variables of the present (k-th) iteration. 

In general, both the prebuckling and buckling variables of equations 
(73), (74), and (75) are represented by Fourier series in the circum- 
ferential coordinate <f>. Since the product of two Fourier series is also 
a Fourier series*, it is seen that each iteration step reduces to an 
ordinary linear statics problem with multiharmonic loading. The solution 
for each component (symmetric and antisymmetric) of each harmonic of the 
effective loading is presented in the preceding discussion of linear 
asymmetric response (SRA 100). Inspection of equations (73), (74), and 
(75) shows that a symmetric prebuckling state results in decoupled 
symmetric and antisymmetric buckling modes, whereas an antisymmetric or 
a general prebuckling state results in a buckling mode with coupled 
symmetric and antisymmetric response (cf. p.17). Even in this case, 
however, the symmetric and antisymmetric components of each harmonic 
for each iteration step are calculated independently (cf. p. 30). 


Inner product .- After each iteration step, the Rayleigh quotient 
[eq. (E-19)] is used to calculate the corresponding eigenvalue estimate. 
For this purpose, it is necessary to be able to compute inner products, 
defined by equation (E-6) . Since prebuckling rotations are neglected, 
equation (E-6) reduces to 


(u;u) = o 0 (1) ’Lii(u,u) - qi(u)' 


u 


(76) 


Evaluation of equation (76) for moderate rotation theory of ring- 
stiffened shells of revolution gives 

(u;u) = /**' {/ (T i ( X X + 06) + T 2 (n + 66) 

< f >=0 

+ T 12 (x^ + ^x) - p(xu + <l<v 

+ (ei + e 2 }w)]rds + £ aT (w w + w w )}d<|> (77) 

j, 9 x x y y 

where u and u are any two kinematically admissible displacement fields, 
the integral over s ranges over the whole meridional length of the shell. 


♦Multiplication of Fourier series is discussed further in the description 
of subroutine MODINT, ref. 8, p. 93. 
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and the summation over r ranges over all attached rings. 

Comparison of equation (76) with the variational eigenvalue 
equation (E-4c) , with terms depending on prebuckling rotations dropped, 
shows that the inner product is, in this case, nothing more than the 
work of the effective loads (for unit X) associated with the displace- 
ment field u (or u) acting through the displacements u (or u) . This 
observation can be made explicit as follows. From equations (1), 
equivalent shell forces are defined as 


Fi = -(r/R 2 )X 1 + r’(X 3 - Lj'/r) 
F 2 = -r’Xi - (r/R 2 )(X 3 - L^/r) 
F 3 = — X 2 + L 3 /R 2 
F 4 = -L 2 


(78) 


and, from equations (21) , equivalent ring forces are defined as 


FLi = -(aF + N ) 
i x y 


FL 2 = _( a F y _ Nx ) 
FL 3 = -(aF^ - N x ) 


(79) 


FL 4 = -aN 


♦ 

Substituting the expressions for u,w,0,i|> from equations (11) and the 
expressions for w x ,Wy from equations (16) into equation (77), then 
performing integrations with respect to <j> by parts, one obtains the 
alternate expression for the inner product 

(u;u) * / 7T [/ Q (F 1 i + F 2 n + F 3 v + F 4 x)rds 

0 s 

+ l (FLiu + FL 2 u + FL 3 u.)]d<f> 


(80) 


where the equivalent forces are given by equations (78) and (79) with 
effective loads given by equations (73), (74), and (75) with X ® 1. 

In equation (80) , each of the displacements and equivalent forces are 
represented by a Fourier series in the circumferential coordinate <j>. 
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Therefore, the integrand and summand are also Fourier series. However, 
because of the <}> -integration only their axisymmetric components contribute 
to the inner product and higher harmonics may be ignored. 


Nonlinear Axisymmetric Response (SRA 200) 


This program solves the nonlinear large-deflection shell equations 
for the case of axisymmetric torsionless loading. The loading is 
assumed to be proportional with the load parameter A, and the first (and 
possibly second) derivatives with respect to A of the response variables 
[i.e. , linear perturbation state(s)] at an input load level Xq, as well 
as the nonlinear response at Ao , are calculated. The numerical solution 
for the nonlinear state is based on a generalization of Newton's method 
for calculating the roots of nonlinear algebraic equations by iteration. 
In addition, for purely mechanical loading, the prebuckling structural 
stiffness Kq (ref. 18) and its derivative at Aq are computed. As 

in reference 18, K is defined as dA/dA, where A is the "work deflection" 
defined such that the area under the A - A curve represents the work of 
the external loading. As shown below, Kq and Kq( 1) are useful in 
calculating the value of a limit load A*, at which the Newton iteration 
does not converge (fig. 5). 


Formulation of equations .- For axisymmetric torsionless loading, shear 
type load and response variables (cf. p. 17 ) are identically zero, and 
the equations associated with them (the axisymmetric torsion equations) 
are dropped from the system of governing equations, thereby reducing their 
differential order from eight to six. Since the general nonlinear 
equations are of the same form as the linear equations plus additional 
effective load terms, the differential equations are obtainable directly 
from equations (22) with n ■ 0. These are 


(rP)' + r[(r/R 2 )X 2 - r'X 3 ) » 0 

(rQ)' - T 2 + r[r'Xx + (r/R 2 )X 3 ] - 0 

(rMi)' + r[r'P - (r/R 2 )Q] - r'M 2 + rL 2 = 0 
, ( 81 ) 
K - r'x - (r/R 2 )ei = 0 

n' + (r/R 2 )x - r'ei = 0 

X* - ki - 0 


The supplemental equations are equations (7a-d), linearized by replacing 
e l> E 2 by ei,e 2 , equations (8a,b), (23a-c), and from equations (23d, e) 
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e 2 


( 82 ) 


= n/r 
k 2 “ r'x/r 

Effective loads, in addition to real proportional (dead) loads XXj, XX 3 , 


X 8 l, X 02 , are from equations ( 2 ), 

L 2 = -Tix (83a) 

Xi - X 3 « Li = 0 (83b) 

plus equation (3a,c,e), (4a,c,e)*, and from equations ( 12 ) 

01 = -X 2 /2 (84a) 

0 2 * 0 (84b) 


Equations (81) and its supplemental equations are a system of six first- 
order differential equations which are of the same form as equation (24), 
where now Y is the six-element column vector (P,Q,Mi ,£ ,ri ,x) and F is a 
nonlinear vector function of Y. Note that all of the nonlinearity is 
exhibited in the two effective load terms L 2 and 81 . 

The effective ring loads due to nonlinear terms [eqs.(20)] are 
identically zero, since the ring rotations are shear type variables. 
Therefore, for axisymmetric torsionless loading, the ring equations and 
associated boundary conditions are linear. Other boundary conditions are 
also assumed to be linear so that equation (25) applies. However, in 
this case [B] and [D] are 3x3 matrices, and the force, displacement , and 
load matrices, {y}, {z}, and (L), are 3x1 matrices. For rings, equations 
(29) -(35) apply with n set to zero and the third row and column of each 
matrix deleted. 

Newton's method .- Following Thurston (ref. 19) a generalization of 
Newton’s method for differential equations is used to reduce the non- 
linear boundary-value problem to a sequence of linear boundary-value 
problems. In this method, the iteration equations are derived by 
assuming that the solution is given by a small correction to an approxi- 
mate solution (Lnitially taken as the linear solution) and linearizing the 


*The effect of the pressure gradient terms inequation (4c) is neglected 
in the program. 
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differential equations with respect to the correction.* Substituting the 
solution after k iterations 



Y 

~ (k— 1) 


+ 6Y 


(85) 


into equation (24), linearizing with respect to 6Y, and using equation (85) 
to eliminate 6Y, gives the iteration equations 



~Y~(k) ~ “ ?Y~(k-l) 


( 86 ) 


where is a matrix, the (i,j) element of which is the derivative of 
i-th component of F with respect to the j-th component of Y, and F and 
F y are both evaluated for Y = Equations (86) are a linear system 

of equations, which, when supplemented by the boundary conditions (25), are 
solved by the Zarghamee method (p. 22). 

For the specialized equations of the previous section, each variable 
can be written as in equation (85). Only the nonlinear terms given by the 
additional loads of equations (83a) and (84a) need be expanded, i.e. 


L *(k) ~ ~ Ti (k) x (k) = _(Tl (k-i) + 6 Ti)(x (k-D + 6x) 

= _Tl (k-l) x (k-l) “ Tl (k-l) 6x " x (k-l) 6Tl 
= _Tl (k-l) x (k) " x (k-l) Tl (k) + Tl (k-l) x (k-l) (87) 

and similarly 

61 (k) ■ -X(k-l) X (k) + (1/2) X(k-l) 2 (88) 


*In this method, each iteration step yields an approximate solution which 
satisfies the boundary conditions exactly, but the differential equations 
only approximately. In reference 20, an alternate form of Newton’s 
method is proposed, whereby the unknowns to be corrected are the initial 
values of Y for each subinterval. In this method, each iteration step 
yields an approximate solution which satisfies the differential equations 
exactly, but the boundary conditions only approximately. 
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The last term on the right-hand side of each of equations (87) and (88) 
are nonhomogeneous terms known from the previous iteration, whereas the 
remaining terms are linear in the variables Tj and x of the present 
iteration. 

It may be observed that in general if * s an exact solution of 

equation (24), then the solution of equation (86) is Y(k) " ?(k-l)» as 
should be in the limit as k °°. However, in this case equations (86) 
represent the variational equations about the equilibrium state Y(k-l)* 
Therefore, if the load Xo is at a limit load X*, a unique solution for 
does not exist (i.e., the system becomes singular). Since a unique 
solution for Y(k)» given an exact solution is a necessary condition 

for convergence, clearly the iteration method diverges at a limit point, 
and in fact, the rate of convergence becomes impractically slow for Xq 
sufficiently close to X*. 

In practice, however, it is not necessary to observe divergence of 
the method in order to estimate X*. This may be done using the stiffness 
Kq and its derivative obtained at Xq < X* (fig. 5). Since at X*, 

Kq = dX/dKo ■ 0, in the vicinity of X*, X may be expanded approximately as 

X » X* + aK 0 2 (89) 


The constant a is evaluated by differentiation of equation (89) with 
respect to Kq to give 


dX/dK 0 



2aK 0 


(90) 


Substitution of equation (90) into equation (89) shows that as Xq 
approaches X*, X* may be computed as 

X* « X 0 - K 0 /2K 0 (1) (91) 

The evaluation of Kq and Kq^ is presented in the description of sub- 
routine STREN in reference 8, p. 106. 

Linear perturbation states .- The calculation of Kq and requires 

the determination of not only the nonlinear response at the given load Xq, 
but also the first and second derivatives of the response (with respect to 
X) at Xq . "f" The differential equations for these states are obtained by 


^The first derivative state is always computed since it is required input 
to the buckling program SRA. 201. 
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differentiation of equation (24) with respect to X to give 

Y (1) ' - F Y Y (1) = 0 (92a) 

Y (2) ' - f y y (2) * Iyy-^ (92b) 

where £yy is a matrix, the (i,j) element of which is the second derivative 
of the i-th component of F with respect to the j-th component of Y, and 
Fy and Fyy are both evaluated at the converged nonlinear solution Y. 

The corresponding boundary conditions are equations (25) with {L} 
replaced by the unit load {L} for the first derivative state and {L} 
replaced by {0} for the second derivative state. Since equations (92) and 
the boundary conditions are linear, they are also solved by the Zarghamee 
method. Comparison of equations (92a, b) with (86) shows that insofar as 
^(k-1) " ?(k) ~ tbe homogeneous forms of these three equations are 
identical. Consequently, the complementary solutions obtained in the last 
Newton iteration may be used in the calculation of y(^) and y(2). 

The nonhomogeneous terms of equations (92a, b) are obtained by 
differentiating with respect to X the dead loads XX^, XX 3 , X0i, X0 2 , the 
live loads given in equations (4a, cl and the effective loads given by 
equations (83a) and (84a). For Y(l), this gives 

Xi = X X + p X 

X 3 = X 3 + p(ei + e 2 ) + £3p/9x + n3p/3y 

L 2 = 0 (93) 

01 - 9l 

0 2 = 0 2 

where C,n,x*ei>e2 are response variables of the nonlinear solution. 
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one obtains 


For Y 


( 2 ) 


Xl 

X 3 

L 2 

8 1 

e 2 


• 2px 


( 1 ) 


2p(ei (1) + e 2 (1) ) + 2£ (1) 3p/3x + 2n (1) 3p/3y 
- 2 Tl (1) X (1) 


"X 

0 


(1)2 


( 94 ) 


Equations (92a) and (92b) are then equivalent to equations (81) and its 
supplemental conditions with nonhomo geneous load terms given by equations 
(93) and (94), respectively. 


Buckling of Axisymmetric Equilibrium States (SRA 201) 

This program calculates the bifurcation buckling inodes of nonlinear 
(or linear) axisymmetric torsionless prebuckling states. The structural 
loading is assumed to have a given spatial distribution, but its magni- 
tude is allowed to vary in proportion to a load parameter X. This leads 
to a set of eigenvalue equations, which are linearized with respect to X 
by expanding the prebuckling state variables in a Taylor series in p = 

X - Xq, and retaining only linear terms in p. 

Geometrically, this method consists of examining the stability of 
fictitious equilibrium states on the tangent to the nonlinear load- 
deformation curve at an assumed load Xq below the critical load. For 
loads near Xq, the corresponding fictitious states are good approximations 
to the neighboring nonlinear states. Consequently, as Xq is increased 
towards the critical load, the fictitious critical loads obtained approxi- 
mate with increasing precision the actual critical load. For each Xq, 
the method of successive approximations is used to obtain the fictitious 
critical load. 

A general form of the eigenvalue equations for bifurcation buckling 
and the method of solution are presented in Appendix E.* In the 
development presented there, no restrictions are placed on the struc- 
tural geometry or loading. In this section the iteration equations 
(E-ll) and the inner product [eq. (E-6) , required to calculate the 
eigenvalue estimate according to the Rayleigh quotient, eq. (E-19)] 
are specialized to ring- stiffened shells of revolution under axisymmet- 
ric torsionless loading. 


*A more specific formulation giving additional details is presented 
in reference 11. 
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Iteration equations .- For numerical purposes, the differential form 
of the variational iteration equations (E-ll) is more convenient to use 
and is derived here. First, the differential form of the variational 
buckling equations (E-l) and eigenvalue equations (E-4) , from which the 
iteration equations follow, are obtained. Applying the usual procedure 
of taking the difference of the governing shell and ring equations 
evaluated for an initial (prebuckling) equilibrium state and an adjacent 
(buckling) equilibrium state, and linearizing the resulting equations 
for the perturbation variables, leads to the basic buckling equations 
[eqs. (E-l)]. The nonlinear shell and ring equations have been presented 
earlier in the form of linear differential equations with nonlinear (and 
live load) terms Isolated as additional effective mechanical and thermal 
loads given by equations (2), (4), (12), (13), (14) and (20). Hence, the 
buckling equations are obtained in the same form as the linear system of 
equations [with dead load terms dropped from the linearized form of eqs. 
(1) , (7), (8) and (21)] with effective additional loads derived from 
equations (2), (4), (12), (13), (14) and (20). However, as mentioned 
previously, the terms of equations (14) are of higher order for moderate 
rotations and are therefore neglected. Furthermore, it is shown in 
reference 11 that it is consistent with this approximation to neglect the 
thermal load of equation (13). For axisymmetric torsionless loading, the 
remaining equations give the following additional effective loads: 


from equations (2) , 

Xi - — t (T i Q + T 2o )6]'/r 
X 2 - -(r'/r) (T 1q + T 2q ) 0 
X 3 = 0 

Ll - T 2q * + x Q Ti2 
L 2 = -<Ti 0 X + X 0 Ti) 

from equations (4) , 

Xi - Xpx 

X2 ■ Apip 

X 3 ■ X[p(e 3 + e 2 ) + £3p/3x 

L| * L»2 * 0 


(95a) 

(95b) 

(95c) 

(95d) 

(95e) 


(96) 

ri3p/9y] 
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from equations (12a) , 


e l = -x 0 x 

02 “ 0 (97) 

®12 = ~X Q ip 

and from equations (20) , 


N = -T , w 

x <Po x 

N = -T . w 

y y (98) 



In equations (95) — (98) , variables with the subscript 0 are prebuckling 
state variables at the load X, and the remaining response variables are 
buckling mode variables. 

Inspection of equations (95) -(98) shows that not only are the equations 
for individual harmonics of the buckling variables decoupled, but also 
that the symmetric and antisymmetric components of each harmonic are 
decoupled. Therefore the buckling equations may be written in terms of 
symmetric or antisymmetric harmonic amplitudes of a single buckling 
harmonic. Equations (22) and their associated supplemental equations 
apply to the amplitudes of symmetric components of the n-th harmonic of 
the buckling mode. The loads for these equations are given by equations 
(3) and (95) — (97) with equation (95a) replaced by the corresponding 
symmetric load amplitude 


Xi » -(n/r)(Ti 0 + T 2() )e 
where from equations (11) 

0 = (r’/r) (v + nq) + (n/R 2 K 


(99) 


( 100 ) 


Boundary conditions for ring boundaries are given by equations (25) , 
(32), (35a), and (35b) with {i t ) = {)l e } = {0} and given by 

equations (30c) and (98) , where from equations (16) the symmetric 
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rotation amplitudes for rings are 


w x = -(nu y + u^)/a 

( 101 ) 

w = nu / a 

y x 


Using the kinematic equation (33), relating shell and ring displace- 
ments, the effective load vector for this case may be written as 


where 


{L} = [<3 [e]{z} 


M 


(T* /a) 


n 2 0 0 

n 2 n 
1 

Symmetric 


0 

0 

0 

0 


(102) 

(103) 


Since the equations for antisymmetric components are the same set of 
equations with n replaced by -n, inspection of these equations shows that 
antisymmetric buckling modes have the same critical value of X as 
symmetric buckling modes, and therefore need not be considered.* 


In order to apply the iteration method, it is necessary to search for 
eigenvalues in a sufficiently small neighborhood of an estimate X = Ao, 
so that in this neighborhood the prebuckling variables Ti , T 2 , Xq , and 

T have a linear dependence on X. Setting 0 ^ 

90 

X = Xq + V (104) 

one has, to first order in 


V x) 

xo.tt) 


T ln + 

0 0 

T 2 + yT2 

0 0 


( 1 ) 


( 1 ) 


xo + PXO 


( 1 ) 


T, + h t a 
<P0 90 


( 1 ) 


(105) 


*Antisymmetric buckling mode shapes are derivable from the corresponding 
symmetric buckling mode shapes simply by changing the sign of the ampli- 
tudes of either the shear or normal type variables. 

Henceforth, variables with the subscript 0 are assumed to be evaluated 
at X ” A q • 
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The eigenvalue problem for P is then obtained by substituting equations 
(104) and (105) into the effective additional load (amplitude) terms 
given by equations (99), (95b-e) , (96), (97), and (102). Each of the 
effective load terms are thus split into two parts, one part independent 
of y, and a second part linear in y. The set of equations so formed 
corresponds to equations (E~4) of Appendix E. 

The iteration equations (E-ll) are obtained by setting y = 1 in the 
second parts of the effective load terms and interpreting the buckling 
mode variables of these parts as being known inputs from the previous 
iteration. Thus the first parts of the effective loads become homogeneous 
terms and the second parts become nonhomo geneous terms for the equivalent 
linear problem of each iteration. 

Inner product .- After each iteration step, the Rayleigh quotient 
[eq. (E-19)] is used to calc u late the corresponding eigenvalue estimate. 
For this purpose, it is necessary to be able to compute inner products, 
defined by equation (E-6).* Evaluation of equation (E-6) for ring- 
stiffened shells of revolution under axisymmetric torsionless loads gives 

(u,o;u,ct) = fll {/[Ti (1) (xx + 00 ) + T 2 (1) (# + 60) 

y u s 0 0 

+ XO^UlX + Tix + Tj. 2 $ + Ti 2 ^P) - PXU - p^v - (p{e + e 2 } 

+ 53p/3x + n3p/3y)w]rds + Y aT, ^ (w w + w w )}d<|> (106) 

r 90 xx y y 

where u and u are any two kinematically admissible displacement fields 
(a and a being corresponding stress states) , the integral over s ranges 
over the whole shell meridian, and the summation over r ranges over all 
attached rings. Since prebuckling variables are axisymmetric, if 
buckling varfables are considered to be amplitudes of . symmetric harmonic 
components, the integral over <p in equation (106) can be replaced by the 
factor it (or 2n in the case of ah axisymmetric buckling mode) . Before 
doing so, however, a more concise form of the inner product is derived 
from equation (106) . 

Comparison of equation (E-6) with equations (E-A) shows that the 
inner product is equivalent to the work of the second part of the 
effective mechanical loads and negative free thermal strains (for unit y) 


*The inner product is also used for orthogonalizing mode estimates with 
respect to lower eigenmodes when calculating nonfundamental eigenvalues. 
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associated with the displacement and stress fields u,a (or u,o) acting 
through the displacement and stress fields u,5 (or u,cr). This observation 
can be made explicit by using equations (78) and (79) and integrating by 
parts with respect to <j>, as was done previously for SRA 101. The result is* 


(u,a;u,cr) - / s ( + F 2 ri + F 3 V + Fpx “ Tl e l “ Ti 2 0i 2 )rds 

+ l (FLiU x + FL 2 u + FL 3 u ) (107) 

r x y 

where the equivalent forces are given by equations (78) and (79) with Li , 
N x * , and replaced by their symmetric components nLi, nN x , and nNy, 
respectively, and effective loads given by equations (99), (95b-e) , (96), 
(97) and (98) , with X replaced by unity and prebuckling quantities 
replaced by derivatives with respect to X at Xq. In equation (107) the 
factor ir (or 2ir for axisymmetric buckling) representing the integration 
with respect to <f> has been dropped. 


Vibrations About Axisymmetric Equilibrium States (SRA 300) 

This program calculates free vibration modes about nonlinear 
axisymmetric torsionless equilibrium states. The eigenvalue equations for 
the square frequency gj 2 of harmonic vibrations about an equilibrium state 
are similar in structure to the eigenvalue equations for the critical load 
increment p for buckling in the vicinity of the same equilibrium state. 

As such they are solved by the same method of successive approximations 
as discussed previously for SRA 201. In this section the iteration 
equations solved and the inner product, used for calculation of the sequence of 
eigenvalue estimates and also for mode orthogonalization when obtaining 
nonfundamental modes, are presented. 

Iteration equations .- The eigenvalue equations for vibrations about 
an equilibrium state are obtained from the eigenvalue equations for 
buckling in the vicinity of the same equilibrium state by replacing the 
incremental loads proportional to p in the buckling problem by the inertial 
loads proportional to to 2 of harmonic vibrations. Based on the thin shell 
assumption that each normal element acts as a rigid body with five degrees 
of freedom, the timewise amplitudes of shell inertial loads due to harmonic 


♦Equation (107) differs in sign (which is immaterial) from the inner 
product as given in reference 11. 
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vibrations of frequency ui are 


*1 

X 2 

X 3 

Ll 

L 2 


o) 2 [m (0) u + m (1) X ] 
m 2 [m^v + 


2 ( 0 ) 

a) z m 'w 


(108) 


-io 2 [m^i|' + m^v] 
ui 2 [m^x + m^u] 


where the mass moments per unit of surface area are given by the integrals 
through the shell wall thickness 


m 


(i) 


/ pz*dz 


for i 0, 1, or 2 


(109) 


Similarly for homogeneous rings, assuming that each cross section acts 
as a rigid element, one obtains 


F - 

tD 2 pAu 


X 

X 


F - 

co 2 pAu 


y 

y 


F <f> “ 

0 J 2 pAu, 

<P 


N = 

w 2 p(I w 

- I w ) 

X 

XX 

xy y 

N = 

0 ) 2 p(l w 

-I w ) 

y 

y y 

xy x' 


w 2 p(I + 
X 

I )w, 

y 4> 


( 110 ) 


As previously discussed for the buckling equations, the eigenvalue 
equations need only be written for the symmetric components of each 
circumferential harmonic, in which case equation (99) replaces equation 
(95a) . The corresponding boundary conditions for ring boundaries are 
given by equations (25) , with [B] and [D] defined by equations (32) , 
and (35a). In addition, the effective load (L) has two components, one 
given by equation (102) and the other obtained from equations (35b) with 
{& t ) = * 0 and {£^( e )> determined by equations (30c) and (110). This 
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second component may be reduced to 


{L} = -(o 2 [y][e]{z} 


(111) 


where 


a 2 A+n 2 I 


n 2 I 


[y] 


£ 

a 


xy 
a 2 A+n 2 I 


nl 


xy 
nl 

x 

a 2 A+I 


0 

0 

0 


Symmetric 


a 2 (I +1 ) 
x y 1 


( 112 ) 


The iteration equations are obtained by setting the eigenvalue 
parameter <o 2 = 1 in equations (108) and ( 111 ) and interpreting u,v,w,x, and 
ip in equations (108) and {z} in equation ( 111 ) as being known inputs from 
the previous iteration. These equations thus give the nonhomo geneous terms 
for the equivalent linear problem of each iteration. 


Inner product .- In analogy with the method presented for buckling 
(see also ref. 10 ), after each iteration for an estimate u(^) of vibration 
mode displacements, an eigenvalue estimate w ^) 2 is obtained from the 
Rayleigh quotient in the form 


w (k ) 2 " (u (k)’ u (k-l) )/(u (k)’ u (k) } (113) 

where (u,u) denotes the inner product of any two kinematically admissible 
displacement fields.* The inner product represents the work of the inertial 
loads associated with the displacement field u (or u) acting through the 
displacement field u ( or u). It may be written in a form similar to 
equation (107). In this case, however, there are no inertial free thermal 
strains and the fourth component of equivalent inertial ring forces 
FL 4 [eq. (79)] is nonzero. Therefore in place of equation (107), one has+ 

(u,u) ■ / g ( F l5 + F 2 n + F 3 v + F 4 x)rds 

+ £ (FL]U x + FL 2 Uy + FI^u^ + FI^w^) (114) 


*The inner product is also used to orthogonalize eigenmode estimates with 
respect to known lower modes in order to force convergence to a higher mode. 

"^Equation (107) differs in sign (which is immaterial) from the definition 
given in reference 10 . 
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where the equivalent forces are given by equations (78) and (79) with Li*, 
N x *, and Ny* replaced by their symmetric components nLi, nN x , and nNy, 
respectively, and effective loads given by equations (108) and (110) with 
u 2 set to unity. 


I CONCLUDING REMARKS 

i 

i 

i 

The governing equations and their method of solution have been pre- 
sented for stress, buckling, and vibration response of branched, stiffened 
shells of revolution tinder axisymmetric and asymmetric loads. In general, 
the numerical solution is reduced to the solution of a sequence of linear 
boundary value problems in ordinary differential equations. These are 
solved by the Zarghamee technique, in which initial conditions for comple- 
mentary and particular solutions (obtained by forward integration) are 
chosen so as to satisfy identically the boundary conditions as the 
corresponding boundaries are passed. This method is more efficient, 
requiring only half as many complementary solutions for open branch 
problems, than the more common method of superposition of complementary 
and particular solutions with arbitrary initial conditions. However, it 
does not eliminate the problem of "long subintervals" associated with 
rapid growth of complementary and particular solutions. 

A further improvement in the method, which eliminates the long 
subinterval problem, as well as providing increased efficiency, is 
currently being evaluated. This new method has been termed the "field 
method," and it eliminates the calculation of complementary and particular 
solutions altogether. The efficiency of the field method may make feasible 
within the scope of the present techniques more advanced problems, such 
as nonlinear response under unsymmetrical loa*ds. 
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APPENDIX A 


The 
given by 


where 


SHELL STIFFNESS COEFFICIENTS 

coefficients X ± , y^ in the shell constitutive equations ( 7 ) are 

Xu - -Ci 2 <0) Al 3 - Ci 2 (1) *14 + C 2 (0> 

X12 " X21 “ “Ci2^ ^13 ” ^12^ ^14 + ^2 

X 1 3 ■ -X31 = Ci2^°^X 33 + C 12 ^X 34 

X14 “ -X41 * Ci2^^X 3 4 + Ci2^ ^44 

X22 “ ”Ci2^X23 “ Cj2^ ^24 + ^2^ 

X23 “ “X 3 2 “ Cj2^^33 + ^12^ ^34 

(A-l) 

X 2 4 * -X 42 * Ci2^X 3 4 + Ci2^ 2 ^44 

X33 “ Ci ( 2) /A 
X34 ■ X43 ■ -Ci ^ /i 

X44 - Ci (0) /A 

un - 4 [G<°>G< 2 > - G< 1 > 2 ]y 2 2 
yi2 m — U2l “ 2[G^ ^ + 2G^ ^/R2]b22 

y 2 2 - [G (0) + 4 G (1) /R 2 + 4 G^ 2 ^/R 2 2 ] -1 

A - Ci<°>Ci< 2 > - Ci 2 (A- 2 ) 
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and for m « 0, 1, or 2 


Cl («) 

* c 8 < m ) + C st < m ) 

c 8 (m) 

- JtEj/U - 

v 1 v 2 )]z m dz 

C 2 (m) 

- J[e 2 /(i - 

viv 2 )]z m dz 

- (m) 

'12 

- /[ V1 E 2 /(1 

- viv 2 )]z in dz 


= Gg( m ) + G st ( m ) 

G s (®) 

= J E 12 z m dz 



In equations (A-3) , the integrals are through the shell wall thickness, 

and C , G ( m ' are stringer contributions given by 
St St 


C st 

(0) 

= NEA/2trr 

c 

(i) 

/-N 

O 
s— ' 

o 

a) 

II 

St 


Z St 

C st 

(2) 

- NEl/2irr + e 2 C 
z 

G 

(0) 

= G ^ =0 

St 


st 

G st 

(2) 

= NGJ/8iTr 


(A-4) 


where e is the normal eccentricity of the centroid of stringer cross 
section! relative to the shell reference surface* 
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APPENDIX B 


SHELLS WITH DOME CLOSURES 


Dome closures are treated by deleting a small spherical cap containing 
the pole and generating appropriate boundary conditions for the artificial 
edge so created.* These boundary conditions, which represent the deleted 
cap to first order in the edge radius, are derived in this appendix. 

As shown in the main body of this report, all of the problems 
treated reduce to the solution of a standard linear statics problem with 
pseudo monoharmonic loads. Furthermore, as noted previously (p. 30 ), the 
solution for antisymmetric load components may be obtained from that for 
equivalent symmetric load problems. Therefore, dome boundary conditions 
need only be written for the case of linear static response under symmetric 
pure harmonic loading. 

Based on the finiteness of strains, one can derive, in a manner 
similar to that of reference 21, the following results valid at a pole 

n 2 £ = 0 
n + nv = 0 

nri + v = 0 (B-l) 

(n 2 -l) x = 0 

k = x 

ei ■ n 

• m 

e 2 = n + nv 
-e 12 = ±nn 
<1 - ±X 
-'P = ±X 

i m x ± n/R 2 


*In evaluation of integrals over the shell meridian, a first-order 
correction is made to account for the contribution over the deleted cap 
(see p. 29 of ref. 13). 


(B-2a) 

(B-2b) 

(B-2c) 

(B-2d) 

(B-2e) 

(B-2f) 
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For the terms with ambiguous signs in equations (B-2) , the lower signs 
apply at a terminal pole, at which r' = -1, and the upper signs apply at 
all other poles, at which r' = 1. The dots above symbols denote 
differentiation with respect to the radius r. Equations (B-l) and (B-2) 
are basic equations from which the dome boundary conditions are derived. 


Zero' th Harmonic (n = 0) 

From equations (B-l) and (B-2) , the following relations hold to 
first order in r at the artificial edge 


n = rn = re 2 

X = rx = ±r<2 


(B-3) 


Using equations (7c, d) and (23c-e) , equations (B-3) become 


±rA 33 Q + rX 34 Mi - (1 + X 13 )n + X 23 x " -r[A 33 0i^ 0 ^ + X34 0 l^] 


(B-4) 


±rX 34 Q + rX 44 Mi - X 14 n + (1 + X 24 ) X = -r[X 34 02 (O) + X 44 02 (1) ] 

Two additional conditions may be derived from the first and third of 
the equilibrium equations (22a) , which for n = 0 reduce to 

(rP)' = r[r'X 3 - (r/R^Xi] 

(B-5) 


(r 2 S) ' = r[(r/R 2 )L 2 - rX 2 ] 

Since ( )' = ( )r', in the vicinity of the pole equations (B-5) may be 
written as 


rP = rX 3 + 0(r 2 ) 
r z S = ±r 2 (Li/R 2 - X 2 ) + 0(r 3 ) 


(B-6) 


where the loads X 3 and L 2 /R 2 - X 2 may be taken as their values at the 
artificial edge. Integration of equations (B-6) between the limits r = 0 
and r = r, and neglecting 0(r 2 ), gives at the edge 


*This is in accordance with the specified description of the reference 
meridian (see pp. 18-19). 
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P - (r/2)X 3 


(B-7) 


iCr/SXL!/^ - X 2 ) 


Equations (B-4) and (B-7) constitute the dome closure boundary 
conditions for n «* 0. They are equivalent to equation (25) with 


[B] 


[D] 


{L} 


1 

0 

0 

0 


0 

±r ^33 

0 

r*34 


0 

0 

1 

0 


0 

±rX 34 

0 

rX 4 4 


0 

0 


0 

0 

0 

(1 + Xi 

3> 

0 

±x 23 

0 

0 


0 

0 

0 

Xi4 


0 ± 

(1 + X 24 ) 

( 

-X 3 /2 





(B-8a) 


(B-8b) 


-r 


X 33©1^°^ + X 34©1^ 

±(X 2 - Li/R 2 )/3 


(B-8c) 


\* 34 0 1 


( 0 ) + 


X4401 


( 1 ) 


First Harmonic (n ■ 1) 

• 

In this case, equations (B-l) yield at the pole £ - rx ■ £ ~ *X 
“ r) + v * 0, so that one may write at a small distance from the pole 

£ - rx ■ 0(r 2 ) (B-9a) 

n + v - r(n + v) + 0(r 2 ) (B-9b) 

Equation (B-9a) yields immediately the first-order condition 

£ - rx - 0 (B-10) 
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However, equation (B-9b) does not yield any new information, since 
substitution of equation (B-2b) for fi + v, and elimination of e 2 through 
equation (23d) results in the identity n + v - p + v. In place of 
equation (B-9b) , equations (B-2a) and (B-2c) can be usftd to form the 
relation 


re* ± rej 2 “ 0(r 2 ) 


(B-ll) 


Using equations (7c, f) and (23c-f) , and neglecting 0(r 2 ), equation (B-ll) 
becomes 


±(*23 ~ V12 ).(5/r - x) - (X13 + X 32 /R 2 )(n + v) ± rX 33 Q ± ry 22 S + rX 34 Mi 

» -r[X 33 0 1 (°> + X 34 0 1 (1) ] (B-12) 

At this point it is desirable to evaluate the term £/r - x* For this 

purpose it is necessary to compute l - rx at the pole. Differentiation 
gives 


5-rx«5-rx- 2x (B-13) 

At the pole r * 0, so that substitution of equation (B-2f) gives 

l - r X - -X ± ri/R 2 (B-14) 


Since, as already noted, £ - rx ■ 5 - rx ■ 0 at the pole, one may write at 
a small distance from the pole 

5 - r X - (r 2 /2)(-x ± n/R 2 ) + 0(r 3 ) (B-15) 

Therefore, to first order in r, the following relation holds 

K/r - x - r(- x ± q/R 2 )/2 « -r(e 12 /R 2 ± <i)/2 (B-16) 


in which equations (B-2c,d) have been used. Substituting equations 
(7d,f) and (23c-f) into equation (B-16) and solving the resulting equation 
for £/r - x gives 


C/r - X - {±(X 14 + X 24 /R 2 ) (n + v)/2 - (r/2) [X 34 Q + p 22 S/R 2 

± X 44 Mx ± X 34 0 1 <°> ± X!*,,©^ 1 )]}/!! + (X 24 - y 12 /R 2 )/2] (B-17) 
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Substitution of equation (B-17) into equation (B-12) to eliminate g/r - x 
gives the desired boundary condition 

±r(Ail33 - A 2 X 3 OQ ± ^22(^1 “ A 2 /R 2 ) S + r(AiX 3 4 - A2X44)M 3 

+ [A-2 (X24 + X24/R2) “ Ai C X 1 3 + X23 /R2) ] (n + v) 

= r[ (A2X34 - A2X33>02^^ + (A2X44 - AiX34)0i^^] (B- 18 ) 


where A* ■ 1 + (X 24 - vi 2 /^- 2 ) 
A2 “ C X 23 “ Mi 2)/2 


(B-19) 


Two additional conditions may be derived from the equilibrium equations 
(22a). For n “ 1, the second and third equations and the first and fourth 
equations may be combined to give 

r’[r(Q - S)J - r[X 2 - r’Xi - (r/R 2 )X 3 ] 


r ' [r(rP + Mi)] - r(r[r'X 3 - (r/R 2 )Xi] - r'L], - L 2 } 
+ (r/R 2 )[r(Q - S) ] 


(B-20) 


In the vicinity 6f the pole, equations (B-20) may be written as 
[r(Q - S)] = r(±X 2 - Xi) + 0(r 2 ) 

(B-21) 

[r(rP + Mi)] = -r(Li ± L 2 ) + 0(r 2 ) 


where the loads Xj , X 2 , Li , L 2 may be taken as their values at the 
artificial edge. Integration of equations (B-21) between the limits r = 0 
and r = r, and neglecting 0(r 2 ) gives at the edge 


Q - S = (r/2) (±X 2 - Xi) 
rP + Mi - -(r/2) (Li ± L 2 ) 


(B-22) 


Equations (B-10) , (B-18) , and (B-22) constitute the dome closure 
boundary conditions for n = 1. They are equivalent to equation (25) with 
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(B-23a) 


[B] 


0 

±r(AiA33-A2A34> 

1 

0 


0 

±ry 22 (Ai-A 2 /R 2 ) 

-1 

0 


0 

r (Aj A34~A 2 A44) 
0 
1 


[D] 


1 0 0 -r 

o a 3 a 3 0 

0 0 0 0 

0 0 0 0 


(B-23b) 


{L} 


l(A 2 A 34 — AiA33 >0 i ^ + (A 2 A 44 ~Ai 134)81 ^ 
(±X 2 -Xi ) /2 
-(L 1 ±L 2 )/2 


(B-23c) 


where 


A3 = A2(A^4 + A24/R2) “ Ai (A 1 3 + A23/R2) (B- 24 ) 

Higher Harmonics (n > 2) 

In this case, equations (B-l) yield at the pole, £ = n=v=x = £ = 0. 
In view of equations (B-2) , the following relations therefore hold to 
first order in r at the artificial edge 

£ = 0 

n = rn = rei 

v = rv = r(ne 2 ± ei 2 )/n 2 (B-25) 

X = rx = ±r<i 

Using equations (7c,d,f) and (23c-f) , equations (B-25) can be put in the 
form of equation (25), where 
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0 

0 

0 

0 

[B] = ± 

0 

±*X 33 

0 

rX 34 


0 

0 

ry 2 2 

0 


0 

±rX 34 

0 

rX 44 


(B-26a) 


[D] 


1 0 

0 -(l+Xi3+n 2 X23/2R2) 

0 ±n(l-yi2/2R2) 

0 -(Xi4+n 2 X24/2R2) 


0 

“*ti( X 1 3+X23/R2) 
0 

-n(Xi4+X24/R2) 


0 

+(l-n 2 /2)X23 

nyi2/2 

+[l+(l-n 2 /2)X 24 ] 


(B-26b) 


{L} ■ -r ■ 


^33©1 (0) + ^34©1 (1) 


^34©1 (0) + *44©1 


(1) 


(B-26c) 
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APPENDIX C 


SUPPORTING LEMMAS FOR ZARGHAMEE METHOD 


In this appendix several basic aspects of the Zarghamee procedure 
are clarified in order to provide some formal justification for this 
method . 


Supplemental Initial Conditions 

As noted in the discussion of the Zarghamee method (p. 24), the 
supplemental initial conditions chosen should insure linear independence 
of the complementary solutions. This question is clarified by the 
following lemma: Any supplemental condition of the form 

[a][U 0 ] + [S][W 0 ] = [I] (C-l) 

where [a] and [$] are 4x4 matrices, insures linear independence of the four 
complementary solutions on open branches. 

The proof is as follows: As is well known, linear independence of a 

set of complementary solutions of a system of ordinary linear differential 
equations implies linear independence everywhere.* Therefore it is 
sufficient to show linear independence at the initial point alone. Linear 
dependence of the complementary solutions at the initial point is equivalent 
to the existence of anon-null 4x1 constant matrix {c} such that 


[boHc} = {0} 
[W 0 ](c} « {0} 


(C-2) 


Postmultiplication of equation (C-l) by (c) shows that equations (C-2) 
imply (c) = {0}, contradicting the requirement of non-null (c). Q.E.D.^ - 


*Since a linear combination of solutions of a system of homogeneous linear 
equations is also a solution, if it vanishes at one point, by the 
uniqueness theorem of initial-value problems, it vanishes everywhere. 

4 > 

This argument may also be used to prove linear independence of the 
additional four complementary solutions on a closed branch. It does not, 
however, prove linear independence of all eight complementary solutions 
on a closed branch. 
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Nonsingularity of [W] and [Z] 


The Zarghamee method involves the inversion of [W] at the final 
points of subintervals on open branches [see eqs.(46) and (49a)], and 
inversion of [W] and [Z] at the final points of subintervals on a closed 
branch [see eqs. (53) and (55)]. In this section it is shown that these 
final point matrices are always nonsingular, so that their inverses exist. 

Let us denote the final value of [W] for a typical subinterval by [Wj ] . 
Consider the columns of the 8x4 matrix of complementary solutions [JJ] as 
the static solutions for a shell whose extent is the subinterval considered 
and which has null surface loads and initial boundary conditions given by 
equation (47a) and the final boundary conditions 

[W] = [W 2 ] (C-3) 


Since equations (C-3) are displacement conditions, it follows from the 
uniqueness theorem of linear elasticity that the corresponding solutions 
are unique. Now assume that [Wi ] is singular, I.e. one of its columns, 
say the i-th column, is a linear combination of the other three.* Since 
equations (47a) and the differential equations are linear and homogeneous, 
and the solutions are unique, it follows that the i-th solution must be 
the same linear combination of the other three solutions. Thus the 
assumption of singular [W i ] contradicts the result of the preceding 
section, in which it is shown that the supplemental initial conditions 
insure the linear independence of complementary solutions. 

The above argument can be applied without change as well to [Z] , for 
which equation (47a) is replaced by equation (54). [For the first sub- 
interval on a closed branch, equations (47a) and (54) are replaced by 
equations (58b).] 


Kinematic Constraints on a Closed Branch 

For the Zarghamee method to be used on a closed branch, it is 
necessary that all eight complementary solutions be linearly independent. 
In this section, the conditions under which eight linearly independent 
complementary solutions exist are derived. 

As noted on pages 23 and 27, the 8x8 matrix of complementary solution 
vectors (k 8) on a closed branch is partitioned into four 


*The possibility that one of its columns is null is excluded since this 
would imply that the corresponding complementary solution is identically 
zero. 
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4x4 .submatrices,* 



(k) 

A necessary and sufficient condition for the Y c to be linearly 
independent is that the determinant of the matrix [Y c ^] be nonzero. As 
noted previously (p. 59), the value of this determinant at any point will 
suffice, and its value A at the initial point of a generic subinterval 
is calculated below. The first subinterval of a closed branch may be 
excluded from further discussion, since in this case it is obvious from 
equations (58a, b) and (C-4) that A = 1. 


For any other subinterval of a closed branch, A is evaluated as 
follows. Assuming that initially U is nonsingular, the first row of the 
right-hand side of equation (C-4) premultiplied by is subtracted 

from the second row to give^ 


u 

V 


u 

V , 

w 

z 


0 

Z-WU V 


“ I U I • I Z-WU _1 V I (C-5) 


In order to evaluate the second factor of this result, note that from 
equation (47a) , 


WIT 1 - -D~^B 


(C-6) 


Substituting equation (C-6) and using equation (54) to eliminate BV, one 
obtains 


Z - Wtf 1 V = Z - D _1 DZ 


D _1 (D - D)Z 


(0-7) 


*For the sake of simplicity, the brackets and braces used for 4x4 and 4x1 
matrices, respectively, are omitted in the remainder of this section. 

^Henceforth in this discussion, it is assumed that U,V,W,Z are evaluated 
at the initial point of the sub interval. 
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Evaluating the difference D-D from equations (48a) , (49a) , and (55) , 
equation (C-7) becomes 


Z - WU _1 V = D^BO^zr 1 - UiWi 1 )Z 


(C-8) 


where, here the subscript 1 denotes values at the end of the preceding 
closed branch subinterval. Substitution of equation (C-8) into equation 
(C-5) gives 

A - |u|-|z|*|B|-|V 1 Zr 1 - UiWi -1 |/|d| (c-9) 

Since from equation (C-6), |w| = -|b|*|u|/|d|, an alternate form of 
equation (C-9) is 

A = - | W| • | Z | • IViZ! -1 - UiWi -1 | (C-10) 


It follows from equation (C-9) that if B is singular, A » 0 and 
therefore linearly independent complementary solutions do not exist in 
the subinterval. Since equation (C-9) is derived without the use of supple- 
mental initial conditions, this conclusion is true regardless of the 
choice of these conditions. Since singular B is equivalent to kinematic 
constraint, for the Zarghamee method kinematic constraints are not allowed 
on a closed branch except at the closure point. 

Equation (C-10) is used to guide the choice of supplemental initial 
conditions on a closed branch. The simplest possible choice consistent 
with the requirement that initially W and Z should be nonsingular is 
used, viz. W = Z = I [cf. eq. (56a)]. 
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APPENDIX D 


CALCULATION OF SHELL STRESSES 


After solution of the shell equations £or the eight basic force and 
displacement variables [fig. 2(a)], SRA 100 (and SRA 200) proceeds to 
calculate all components of the three-dimensional stress tensor except 
the transverse normal stress cr z , which is negligible for a thin shell. 

The equations used for the calculation of these stresses are presented 
in this appendix. 

According to the Love-Kirchhoff thin shell hypothesis, neglecting 
terms of order z/R relative to unity,* the primary (i.e. , in-surface) 
shell stresses for an orthotropic wall are given by 

°s = tEi / (1 - v 1 v 2 )][ei + v 2 e 2 + Z (<1 + v 2 k 2 ) - (0! + v 2 0 2 )] 
a 4 , ■ [E 2 /(l - Viv 2 )][e 2 + Vi£! + z(k 2 + Vitq) - (0 2 + v^)] (D-l) 

a s4> “ E 1 2 ( G 1 2 + 2zk 1 2 ) 

In order to compute these stress components from eqs. (D-l), all that is 
needed are the stretching (ei,e2» e 12) and bending (ki,K2> k 12) strains of 
the reference surface. After solution of the differential equations for 
the basic force (y) and displacement (z) vectors (see pp. 22-30), these 
strains are given by eqs. (23d-f) and (7c,d,f) with the aid of eqs. (23c) 
and (8) . 

The transverse shear and normal stresses may then be obtained by 
integration through the shell thickness of the three-dimensional equa- 
tions of equilibrium. In terms of symmetrical stress amplitudes (table 
I) these equations are, assuming no body forces and neglecting terms of 
order z/R relative to unity, 


9o sz /3z = -3o g /9s - (n/r )a a< p + (r'/r)(a (() - o s ) (D-2a) 

3o ( j )Z /3z = -gog^/gs + (n/r)^ - 2(r , /r)o S( j ) (D-2b) 

9a z /3z = -(l/r)3(ro sz )/3s - (n/r )a^ 2 + a g / Ri + a^/ R 2 (D-2c) 

where Rj is the meridional radius of curvature. Inspection of eqs. (D-2) 
shows that for thin shells the transverse stresses o sz ,o ( p z ,o z can be 


*Here, R represents either principal radius of curvature (Rj or R 2 ) of 
the reference surface. 
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significant (i.e., >> z/R times a primary stress) only in a boundary-layer 
zone, where s-differentiation is tantamount to multiplication by a factor 
much larger than 1/R, or if n/r >> 1/R. On the other hand, even in a 
boundary-layer zone, for the thin shell hypothesis to remain valid it is 
necessary that c sz ,c<j, z >c z << c s >°<|>» 0 s<j>* Accepting this as being the case, 
it is clear from eqs. (D-2) that the transverse normal stress a z is 
negligible, even in boundary-layer zones and hence need not be computed. 

In practice, eqs. (D-2a,b) are integrated with respect to z starting 
at the shell wall inner face. At this face o sz and o^ z are given as -Xj 
and -X 2 by eqs. (3a, b) for an attached elastic foundation or are simply 
zero if no foundation exists. In evaluating the s-derivatives on the 
right-hand sides of eqs. (D-2a,b) the differential equations (22) are 
used to obtain derivatives of response variables, whereas derivatives 
of wall properties and thermal loads are neglected. 
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APPENDIX E 


GENERAL BUCKLING EQUATIONS 


In this appendix the treatment of the buckling equations (based on 
moderate rotation theory) without restriction of structural geometry or 
loading is presented. This development may therefore be taken as the 
basis for both the bifurcation analysis of linearized asymmetric 
equilibrium states (SRA 101) and the bifurcation analysis of nonlinear 
axisymmetric equilibrium states (SRA 201) . Although set in a more 
general context, this development follows along the same lines as that 
of reference 11. 


Eigenvalue Equations 

Buckling equations without restriction on structural geometry or 
loading distribution have been presented in reference 18. As presented 
there, the variational form of these equations for the buckling mode 
variables u, a , e in terms of the prebuckling equilibrium state 
variables uq(X), oq(X) at a load factor X (for a proportional loading) is 

e » Li(u) + L n (u 0 ,u) 

a = H(e) (E-l) 

o*6e 0 + ao*Ln(u,6u) - Aqi(u)*6u = 0 


where 


6e 0 = L x (6u) + Ln(u 0 ,6u) (E-2) 

In order to apply the iteration method for bifurcation problems, it is 
necessary to search for eigenvalues in a sufficiently small neighborhood 
of a load estimate X = Xq, so that in this neighborhood the prebuckling 
state variables uq(X), Oq(A) have a linear dependence on X. Setting 
X - Xq + y, one has, to first order in y,^" 

uo(A) = uq + yuo^ 

m (E ~ 3) 

oq(X) = ctq + yoo'' ' 

Substitution of equations (E-3) into equations (E-l) and (E-2) yields the 


tHenceforth, variables with the subscript 0 are assumed to be evaluated 
at X = Xq. 
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linear eigenvalue problem with eigenvalues p^, i = 1,2, ... (arranged in the 
order of increasing absolute value) 

e - LjCu) + LuCuqju) + pLuCuq^ 1 ^ ,u) (E-4a) 

o ■ H(e) (E-4b) 

a*6eo + ao*Ln(u,6u) - Xoqi(u)*6u 

+ y[a*Lii(uo^^,6u) + Oq ^ *Li i (u,6u) - q 1 (u)*6u] - 0 (E-4c) 

Orthogonality of eigenfunctions is obtained by setting 6u * u. in 
equations (E-4c) evaluated for the i-th eigenvalue, interchanging 1 and 
j, and subtracting the resulting two equations. With the aid of 
equations (E-2) , (E-4a) and (E-4b) , this gives the result 

(y i - Pj)(u 1 ,o i ;Uj,aj) = 0 (E-5) 

where the inner (scalar) product is defined by 

(u i ,a 1 ;uj,0j) 5 Oo^'LnO^.Uj) + o^Ln (u 0 (1 \uj) 

+ Oj-Lutuo^ 1 ^ ,u ± ) - qi (u ± ) • Uj (E-6) 

It is noted that the second and third terms on the right-hand side of 
equation (E-6) represent prebuckling deformations (i.e., rotations), and 

the last term represents conservative live loads (i.e., a normal pressure 

field). Also, the inner product is seen to be a symmetric functional, 
i.e., (u^c^jUjjOj) - (Uj ,0^ ^jC^). From equation (E-5) it follows that 

eigenfunctions corresponding to distinct eigenvalues p. ,p. are orthogonal, 
i.e. 3 


(u i ,a i ;u^,Oj) =0 , i ^ j (E-7) 

Similarly, evaluation of equation (E-4c) for the i-th eigenvalue with 
6u * u^ gives the expression for the eigenvalues in terms of the 
eigenfunctions, viz. 

vt ± “ ~l a o *1*2 (u ± ) + a i* e i " ^oqi(u i )*u i ]/[a 0 ( 1 ^*L 2 (u i ) - qi (u ± > *u ± ) (E-8) 


*For X q ■ X c , the critical load, this definition of the inner product is 
equal to the functional f(1)(u^,Uj) defined in reference 18. 


66 


The right-hand side of equation (E-8) (with the index i deleted) defines 
the Rayleigh quotient, which may be used to obtain an estimate y for y^ 
if estimates u, a, e for Ui, 02 , ei are available. However, in applying 
equation (E-8) in this way, care must be taken that u, e and y satisfy the 
kinematic equation (E-4a) , in addition to a, e satisfying equation (E-4b) . 
Such kinematically admissible functions admit the eigenfunction expansions 


u " l Vi 

i-1 


0 - I Vi 

i=l 1 1 


(E-9) 


Forming from equations (E-9) the inner product (u,o;u. ,a.) shows that, in 
view of equation (E-7), ^ ^ 

k ± = (u,o;u i ,a i )/(u 1 ,o i ;u i ,o i ) (E-10) 


Iterative Solution of Eigenvalue Equations 

The iteration method consists of successive solution of the following 
modified form of equations (E-4) 


e (k) = Ll(u (k)^ + L H (u 0» u ( k )) + L ll< u 0 (1) » u ( k _ 1 )> (E-lla) 

a (k) " H < e (k) ) (E ' llb) 

a (k)‘ 6e ° + a 0 * L 1 1 ( u ( k ) » 6u ) “*0<!l< u ( k )) #fiu + CT( k _i) * L 1 1 *<Su> 

+ o 0 ^ 1 ^Li 1 (u^ k _ 1 ^,6u) - q 1 (u^ k _ 1 ^)*6u - 0 (E-llc) 


where the subscript (k) denotes the solution after k iterations. The 
estimates u( k \ , a ( k )» e (k) cannot bs used directly in the Rayleigh 
quotient to obtain an estimate U( k )* nor d° they admit eigenfunction 
expansions, since in general U( k \ , e^s, and y,,.. would not satisfy the 
kinematic equations (E-4a) . It is therefore necessary to relate 


.(H) ’ e (k) 


to modified variables a 


and y^ k ^ are kinematically admissib 




such that u 


(k) * °(k)» e (k) 


e (k) " Ll(u (k) ) + L H ( u 0> u ( k )) + v (k) Lll(u ° (1 ^ » u (k)^ (E-12a) 

*(k) = H(e (k) ) (E-12b) 
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Taking the difference of equations (E-12a,b) and (E-lla,b) gives 


" V> - H( '( k ) - e (k)> 

£ (k) • E (k) ■ L n ( “» C1> ’ u (k) u (k) - "(k-l) ) 


A 

from which it is observed that the difference - 0( k \ is of the order 

of prebuckling rotations, the square of which is negligible for our 
moderate rotation equations. Forming the work term 0(k)* e (k) (required 
in the Rayleigh quotient) from equations (E-13) , and neglecting terms of 
the order of the square of prebuckling rotations, gives 

“<k) - '(w ■ °(k) -E (k) + 2a (k)’ L > i(u ° (1) -‘ , <k)"(k) • °(k-i) > <E_U> 

Setting 6u = u^ in equation (E-llc) , 6u = U( k ) in equation (E-4c) 
evaluated for the i-th eigenvalue, and subtracting the results gives the 
relation 


(u i ,a i ;u (k) , °(k) ) “ ( V a i ;u (k-l)’ a (k-l) )/y i 

From equation (E-6) it follows that 

( u i»°i > u (k) ,a (k)> “ (u i’ a i ;u (k) ,a (k) ) 

■ <°( k ) _ ° (k) ) , L 11 ( u 0 <1) .« 1 ) 


(E-15) 


(E-16) 


Since as observed previously a/. \ - a,, is of the order of prebuckling 
rotations, this difference of inner products is of the order of the 
square of prebuckling rotations and therefore negligible. Hence, it 
follows from equations (E-10) and (E~15) that the expansion coefficients 
f° r t ^ ie functions U( k ) , cr^) satisfy the relation 


A., = A. n , /y . 
lk 1 ,k-l 1 


(E-17) 


from which it follows that as k 00 


U (k) - 



proving convergence of the iteration method. 


(E-18) 
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Substituting U(k) , ° (v) » e (k) u i 5 ^ a i» e i dn fc h e Rayleigh quotient 
given in equation (E-8) , eliminating a (k)‘ e (k) through use of equation 
(E-14), eliminating a (k)* e (k) through use or equation (e-11c) with 
<$u = u (k) * and solving for the corresponding eigenvalue estimate P( k ) 
gives the result 


00 


= (u 


(k) 


,0 (k) ;u (k-l) ,<7 (k-l) )/(u (k) ,a (k) ;u (k) > °(k) ) 


(E-19) 


Since, as already noted, the inner product is insensitive to the difference 
CT {k) “ °(k)’ to the accurac y of our moderate rotation equations, P( k ) is 
also given by equation (E-19) with o^) replaced by °^k)* this form can 
be applied the completeness relation 


(u (k) , °(k) ;u (k-l) , 0 (k-l) ) J ( V a i ; V a i )A ik A i,k-l 

I s *! 

giving in view of equation (E-17) the result 

QO 00 

"oo "Jji'VVVVV'^VVvVV 


(E-20) 


(E-21) 


Equation (E-21) shows that, by use of the Rayleigh quotient, P(k) 
converges to iq at a much faster rate than U( k ) > 0(k) converge to uj, oj. 


Finally, since p( k ) "*■ Mi» and from equation (E-17) -*• u (k-i') > 

it follows from equations (E-13) that cJ(k) -*■ a/ k \. Hence, in spite ot J 
its inadmissibility 


(k) 


oi 


(E-22) 


Thus, in the practical application of the iteration procedure, at no time 
is it necessary to work with the modified variables a^) ° r e (k)* 
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TABLE I 


CIRCUMFERENTIAL VARIATION OF HARMONIC VARIABLES 


Symmetric 


Antisymmetric 


cos n<J> 


sin n<J> 


cos n<j> 


ariables 



Xi,X 3 

l 2 

01 

01,02 


F X ,Fy 



N x ,Ny 


N x ,Ny 


F x » F y 

Na 


Ti,T 2 

Mi,M 2 


Ti,T 2 

Mi,M 2 


Ring 

Response 

Variables 


e l » e 2 
K 1 » K 2 

a s» a 4>» 0 z» a sz 


T * 

Mx.My 

u x ,u y 


ei.^2 
K 1 » K 2 


Mx.My 

u x> u y 


K x ,K y 


K x» K y 



















AXIS OF REVOLUTION 




E: 


The reference surface may be chosen as any convenient 
continuous surface within or near the shell wall. 


FIGURE 1. HYPOTHETICAL BRANCHED SHELL PROFILE 

(WITH FIVE BRANCH POINTS AND ONE CLOSED BRANCH) 




AXIS OF REVOLUTION 



FIGURE 2. POSITIVE DIRECTIONS FOR SHELL DISPLACEMENTS, FORCES, AND LOADS 
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NOTE: Coordinates s,<|>,z form a right-handed system. The radius of curvature 

r 2 is reckoned positive only if the local positive z-direction points 
away from the axis of revolution. 



RING CENTRO IDAL AXIS 



FIGURE 3. POSITIVE DIRECTIONS FOR RING VARIABLES 




BRANCH POINT 



NOTE: Arrows indicate direction of increasing s. 

Subintervals ij, ± 2 j i3, and 14 (here J = 4 ) 
enter the branch point, whereas subinterval 
i4+l exits the branch point. The branch 
point is denoted as ij,l on entering sub- 
intervals and ij+1,0 on the exiting subinterval. 


FIGURE 4 . 


HYPOTHETICAL BRANCH POINT BOUNDARY 
(WITH FOUR ENTERING SUBINTERVALS) 


A 


A 



A = 0 to A = Aq 


(a) LOAD VS. WORK DEFLECTION 



LINEAR STIFFNESS 


(b) LOAD VS. STIFFNESS 
K = dA/dA 


FIGURE 5. TYPICAL PREBUCKLING LOADING DIAGRAMS 


NASA-Langley, 1973 32 CR-2085 
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